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

Merge pull request #3996 from lafourcadep/snann_slcsa

Browse Source 
Compute sna/atom on fixed number of neighbors and compute slcsa/atom (Supervised Learning Crystal Structure Analysis tool)
This commit is contained in:
Axel Kohlmeyer
2023-12-13 14:36:49 -05:00
committed by GitHub
parent 0c4a1cb21d 754041ee36
commit 89c5ee7079
21 changed files with 2559 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
doc/src/Bibliography.rst
View File

@ -562,6 +562,9 @@ Bibliography
**(Kumar)**
Kumar and Skinner, J. Phys. Chem. B, 112, 8311 (2008)
**(Lafourcade)**
Lafourcade, Maillet, Denoual, Duval, Allera, Goryaeva, and Marinica, `Comp. Mat. Science, 230, 112534 (2023) <https://doi.org/10.1016/j.commatsci.2023.112534>`_
**(Lamoureux and Roux)**
G.\ Lamoureux, B. Roux, J. Chem. Phys 119, 3025 (2003)

1
doc/src/Commands_compute.rst
View File

@ -123,6 +123,7 @@ KOKKOS, o = OPENMP, t = OPT.
* :doc:`reduce/region <compute_reduce>`
* :doc:`rigid/local <compute_rigid_local>`
* :doc:`saed <compute_saed>`
* :doc:`slcsa/atom <compute_slcsa_atom>`
* :doc:`slice <compute_slice>`
* :doc:`smd/contact/radius <compute_smd_contact_radius>`
* :doc:`smd/damage <compute_smd_damage>`

1
doc/src/compute.rst
View File

@ -287,6 +287,7 @@ The individual style names on the :doc:`Commands compute <Commands_compute>` pag
* :doc:`reduce/region <compute_reduce>` - same as compute reduce, within a region
* :doc:`rigid/local <compute_rigid_local>` - extract rigid body attributes
* :doc:`saed <compute_saed>` - electron diffraction intensity on a mesh of reciprocal lattice nodes
* :doc:`slcsa/atom <compute_slcsa_atom>` - perform Supervised Learning Crystal Structure Analysis (SL-CSA)
* :doc:`slice <compute_slice>` - extract values from global vector or array
* :doc:`smd/contact/radius <compute_smd_contact_radius>` - contact radius for Smooth Mach Dynamics
* :doc:`smd/damage <compute_smd_damage>` - damage status of SPH particles in Smooth Mach Dynamics

162
doc/src/compute_slcsa_atom.rst Normal file
View File

@ -0,0 +1,162 @@
.. index:: compute slcsa/atom
compute slcsa/atom command
============================
Syntax
""""""
.. code-block:: LAMMPS
compute ID group-ID slcsa/atom twojmax nclasses db_mean_descriptor_file lda_file lr_decision_file lr_bias_file maha_file value
* ID, group-ID are documented in :doc:`compute <compute>` command
* slcsa/atom = style name of this compute command
* twojmax = band limit for bispectrum components (non-negative integer)
* nclasses = number of crystal structures used in the database for the classifier SL-CSA
* db_mean_descriptor_file = file name of file containing the database mean descriptor
* lda_file = file name of file containing the linear discriminant analysis matrix for dimension reduction
* lr_decision_file = file name of file containing the scaling matrix for logistic regression classification
* lr_bias_file = file name of file containing the bias vector for logistic regression classification
* maha_file = file name of file containing for each crystal structure: the Mahalanobis distance threshold for sanity check purposes, the average reduced descriptor and the inverse of the corresponding covariance matrix
* c_ID[*] = compute ID of previously required *compute sna/atom* command
Examples
""""""""
.. code-block:: LAMMPS
compute b1 all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1 delta 0.3
compute b2 all slcsa/atom 8 4 mean_descriptors.dat lda_scalings.dat lr_decision.dat lr_bias.dat maha_thresholds.dat c_b1[*]
Description
"""""""""""
.. versionadded:: TBD
Define a computation that performs the Supervised Learning Crystal
Structure Analysis (SL-CSA) from :ref:`(Lafourcade) <Lafourcade2023_1>`
for each atom in the group. The SL-CSA tool takes as an input a per-atom
descriptor (bispectrum) that is computed through the *compute sna/atom*
command and then proceeds to a dimension reduction step followed by a
logistic regression in order to assign a probable crystal structure to
each atom in the group. The SL-CSA tool is pre-trained on a database
containing :math:`C` distinct crystal structures from which a crystal
structure classifier is derived and a tutorial to build such a tool is
available at `SL-CSA <https://github.com/lafourcadep/SL-CSA>`_.
The first step of the SL-CSA tool consists in performing a dimension
reduction of the per-atom descriptor :math:`\mathbf{B}^i \in
\mathbb{R}^{D}` through the Linear Discriminant Analysis (LDA) method,
leading to a new projected descriptor
:math:`\mathbf{x}^i=\mathrm{P}_\mathrm{LDA}(\mathbf{B}^i):\mathbb{R}^D
\rightarrow \mathbb{R}^{d=C-1}`:
.. math::
\mathbf{x}^i = \mathbf{C}^T_\mathrm{LDA} \cdot (\mathbf{B}^i - \mu^\mathbf{B}_\mathrm{db})
where :math:`\mathbf{C}^T_\mathrm{LDA} \in \mathbb{R}^{D \times d}` is
the reduction coefficients matrix of the LDA model read in file
*lda_file*, :math:`\mathbf{B}^i \in \mathbb{R}^{D}` is the bispectrum of
atom :math:`i` and :math:`\mu^\mathbf{B}_\mathrm{db} \in \mathbb{R}^{D}`
is the average descriptor of the entire database. The latter is computed
from the average descriptors of each crystal structure read from the
file *mean_descriptors_file*.
The new projected descriptor with dimension :math:`d=C-1` allows for a
good separation of different crystal structures fingerprints in the
latent space.
Once the dimension reduction step is performed by means of LDA, the new
descriptor :math:`\mathbf{x}^i \in \mathbb{R}^{d=C-1}` is taken as an
input for performing a multinomial logistic regression (LR) which
provides a score vector
:math:`\mathbf{s}^i=\mathrm{P}_\mathrm{LR}(\mathbf{x}^i):\mathbb{R}^d
\rightarrow \mathbb{R}^C` defined as:
.. math::
\mathbf{s}^i = \mathbf{b}_\mathrm{LR} + \mathbf{D}_\mathrm{LR} \cdot {\mathbf{x}^i}^T
with :math:`\mathbf{b}_\mathrm{LR} \in \mathbb{R}^C` and
:math:`\mathbf{D}_\mathrm{LR} \in \mathbb{R}^{C \times d}` the bias
vector and decision matrix of the LR model after training both read in
files *lr_fil1* and *lr_file2* respectively.
Finally, a probability vector
:math:`\mathbf{p}^i=\mathrm{P}_\mathrm{LR}(\mathbf{x}^i):\mathbb{R}^d
\rightarrow \mathbb{R}^C` is defined as:
.. math::
\mathbf{p}^i = \frac{\mathrm{exp}(\mathbf{s}^i)}{\sum\limits_{j} \mathrm{exp}(s^i_j) }
from which the crystal structure assigned to each atom with descriptor
:math:`\mathbf{B}^i` and projected descriptor :math:`\mathbf{x}^i` is
computed as the *argmax* of the probability vector
:math:`\mathbf{p}^i`. Since the logistic regression step systematically
attributes a crystal structure to each atom, a sanity check is needed to
avoid misclassification. To this end, a per-atom Mahalanobis distance to
each crystal structure *CS* present in the database is computed:
.. math::
d_\mathrm{Mahalanobis}^{i \rightarrow \mathrm{CS}} = \sqrt{(\mathbf{x}^i - \mathbf{\mu}^\mathbf{x}_\mathrm{CS})^\mathrm{T} \cdot \mathbf{\Sigma}^{-1}_\mathrm{CS} \cdot (\mathbf{x}^i - \mathbf{\mu}^\mathbf{x}_\mathrm{CS}) }
where :math:`\mathbf{\mu}^\mathbf{x}_\mathrm{CS} \in \mathbb{R}^{d}` is
the average projected descriptor of crystal structure *CS* in the
database and where :math:`\mathbf{\Sigma}_\mathrm{CS} \in \mathbb{R}^{d
\times d}` is the corresponding covariance matrix. Finally, if the
Mahalanobis distance to crystal structure *CS* for atom *i* is greater
than the pre-determined threshold, no crystal structure is assigned to
atom *i*. The Mahalanobis distance thresholds are read in file
*maha_file* while the covariance matrices are read in file
*covmat_file*.
The `SL-CSA <https://github.com/lafourcadep/SL-CSA>`_ framework provides
an automatic computation of the different matrices and thresholds
required for a proper classification and writes down all the required
files for calling the *compute slcsa/atom* command.
The *compute slcsa/atom* command requires that the :doc:`compute
sna/atom <compute_sna_atom>` command is called before as it takes the
resulting per-atom bispectrum as an input. In addition, it is crucial
that the value *twojmax* is set to the same value of the value *twojmax*
used in the *compute sna/atom* command, as well as that the value
*nclasses* is set to the number of crystal structures used in the
database to train the SL-CSA tool.
Output info
"""""""""""
By default, this compute computes the Mahalanobis distances to the
different crystal structures present in the database in addition to
assigning a crystal structure for each atom as a per-atom vector, which
can be accessed by any command that uses per-atom values from a compute
as input. See the :doc:`Howto output <Howto_output>` page for an
overview of LAMMPS output options.
Restrictions
""""""""""""
This compute is part of the EXTRA-COMPUTE package. It is only enabled
if LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
Related commands
""""""""""""""""
:doc:`compute sna/atom <compute_sna_atom>`
Default
"""""""
none
----------
.. _Lafourcade2023_1:
**(Lafourcade)** Lafourcade, Maillet, Denoual, Duval, Allera, Goryaeva, and Marinica,
`Comp. Mat. Science, 230, 112534 (2023) <https://doi.org/10.1016/j.commatsci.2023.112534>`_

39
doc/src/compute_sna_atom.rst
View File

@ -45,7 +45,7 @@ Syntax
* w_1, w_2,... = list of neighbor weights, one for each type
* nx, ny, nz = number of grid points in x, y, and z directions (positive integer)
* zero or more keyword/value pairs may be appended
* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* or *dgradflag*
* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* or *dgradflag* or *nnn* or *wmode* or *delta*
.. parsed-literal::
@ -82,6 +82,16 @@ Syntax
*0* = descriptor gradients are summed over atoms of each type
*1* = descriptor gradients are listed separately for each atom pair
* additional keyword = *nnn* or *wmode* or *delta*
.. parsed-literal::
*nnn* value = number of considered nearest neighbors to compute the bispectrum over a target specific number of neighbors (only implemented for compute sna/atom)
*wmode* value = weight function for finding optimal cutoff to match the target number of neighbors (required if nnn used, only implemented for compute sna/atom)
*0* = heavyside weight function
*1* = hyperbolic tangent weight function
*delta* value = transition interval centered at cutoff distance for hyperbolic tangent weight function (ignored if wmode=0, required if wmode=1, only implemented for compute sna/atom)
Examples
""""""""
@ -94,6 +104,7 @@ Examples
compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 chem 2 0 1
compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1 sinner 1.35 1.6 dinner 0.25 0.3
compute bgrid all sna/grid/local 200 200 200 1.4 0.95 6 2.0 1.0
compute bnnn all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1 delta 0.2
Description
"""""""""""
@ -433,6 +444,25 @@ requires that *bikflag=1*.
The rerun script can use a :doc:`special_bonds <special_bonds>`
command that includes all pairs in the neighbor list.
The keyword *nnn* allows for the calculation of the bispectrum over a
specific target number of neighbors. This option is only implemented for
the compute *sna/atom*\ . An optimal cutoff radius for defining the
neighborhood of the central atom is calculated by means of a dichotomy
algorithm. This iterative process allows to assign weights to
neighboring atoms in order to match the total sum of weights with the
target number of neighbors. Depending on the radial weight function
used in that process, the cutoff radius can fluctuate a lot in the
presence of thermal noise. Therefore, in addition to the *nnn* keyword,
the keyword *wmode* allows to choose whether a Heaviside (*wmode* = 0)
function or a Hyperbolic tangent function (*wmode* = 1) should be used.
If the Heaviside function is used, the cutoff radius exactly matches the
distance between the central atom an its *nnn*'th neighbor. However, in
the case of the hyperbolic tangent function, the dichotomy algorithm
allows to span the weights over a distance *delta* in order to reduce
fluctuations in the resulting local atomic environment fingerprint. The
detailed formalism is given in the paper by Lafourcade et
al. :ref:`(Lafourcade) <Lafourcade2023_2>`.
----------
Output info
@ -585,6 +615,7 @@ Related commands
""""""""""""""""
:doc:`pair_style snap <pair_snap>`
:doc:`compute slcsa/atom <compute_slcsa_atom>`
Default
"""""""
@ -592,6 +623,7 @@ Default
The optional keyword defaults are *rmin0* = 0,
*switchflag* = 1, *bzeroflag* = 1, *quadraticflag* = 0,
*bnormflag* = 0, *wselfallflag* = 0, *switchinnerflag* = 0,
*nnn* = -1, *wmode* = 0, *delta* = 1.e-3
----------
@ -623,3 +655,8 @@ of Angular Momentum, World Scientific, Singapore (1987).
.. _Ellis2021:
**(Ellis)** Ellis, Fiedler, Popoola, Modine, Stephens, Thompson, Cangi, Rajamanickam, Phys Rev B, 104, 035120, (2021)
.. _Lafourcade2023_2:
**(Lafourcade)** Lafourcade, Maillet, Denoual, Duval, Allera, Goryaeva, and Marinica,
`Comp. Mat. Science, 230, 112534 (2023) <https://doi.org/10.1016/j.commatsci.2023.112534>`_

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

@ -79,6 +79,7 @@ Alessandro
Alexey
ali
aliceblue
Allera
Allinger
allocatable
allocator
@ -719,6 +720,7 @@ dem
Dendrimer
dendritic
Denniston
Denoual
dephase
dephasing
dequidt
@ -860,6 +862,7 @@ Dunweg
Dupend
Dupont
dUs
Duval
dV
dvector
dVx
@ -1299,6 +1302,7 @@ Gonzalez-Melchor
googlemail
googletest
Gordan
Goryaeva
Goudeau
GPa
GPL
@ -1385,6 +1389,7 @@ hcp
hdnnp
HDNNP
Hearn
Heaviside
heatconduction
heatflow
Hebbeker
@ -1849,6 +1854,7 @@ lbl
LBtype
lcbop
ld
lda
ldfftw
ldg
lebedeva
@ -1970,6 +1976,7 @@ lossy
Lozovik
lps
lpsapi
lr
lrt
lsfftw
lt
@ -2012,7 +2019,10 @@ magelec
Maginn
magneton
magnetons
maha
Mahalanobis
Mahoney
Maillet
mainboard
mainboards
makefile
@ -2191,6 +2201,7 @@ mintcream
Mintmire
Miron
mis
misclassification
Mises
Mishin
Mishra
@ -2299,6 +2310,7 @@ multicomponent
multicore
multielectron
multinode
multinomial
multiphysics
Multipole
multiscale
@ -2390,6 +2402,7 @@ Nbtypes
Nbytes
nc
Nc
nclasses
nchunk
Nchunk
ncoeff
@ -3352,6 +3365,7 @@ Skylake
slateblue
slategray
slater
slcsa
Slepoy
Sliozberg
sLL

1
examples/PACKAGES/sna_nnn_slcsa/Zr_mm.eam.fs Symbolic link
View File

@ -0,0 +1 @@
../../../potentials/Zr_mm.eam.fs

879
examples/PACKAGES/sna_nnn_slcsa/data.zr_cell Normal file
View File

@ -0,0 +1,879 @@
# Hcp Zr with box vectors H1=[2-1-10], H2=[-12-10], H3=[0001].
864 atoms
1 atom types
0.000000000000 19.374000000000 xlo xhi
0.000000000000 33.556752345839 ylo yhi
0.000000000000 30.846000000000 zlo zhi
Masses
1 91.22400000 # Zr
Atoms # atomic
1 1 0.000000000000 1.864264019213 2.570500000000
2 1 0.000000000000 0.000000000000 0.000000000000
3 1 1.614500000000 4.660660048033 2.570500000000
4 1 1.614500000000 2.796396028820 0.000000000000
5 1 3.229000000000 1.864264019213 2.570500000000
6 1 3.229000000000 0.000000000000 0.000000000000
7 1 4.843500000000 4.660660048033 2.570500000000
8 1 4.843500000000 2.796396028820 0.000000000000
9 1 6.458000000000 1.864264019213 2.570500000000
10 1 6.458000000000 0.000000000000 0.000000000000
11 1 8.072500000000 4.660660048033 2.570500000000
12 1 8.072500000000 2.796396028820 0.000000000000
13 1 9.687000000000 1.864264019213 2.570500000000
14 1 9.687000000000 0.000000000000 0.000000000000
15 1 11.301500000000 4.660660048033 2.570500000000
16 1 11.301500000000 2.796396028820 0.000000000000
17 1 12.916000000000 1.864264019213 2.570500000000
18 1 12.916000000000 0.000000000000 0.000000000000
19 1 14.530500000000 4.660660048033 2.570500000000
20 1 14.530500000000 2.796396028820 0.000000000000
21 1 16.145000000000 1.864264019213 2.570500000000
22 1 16.145000000000 0.000000000000 0.000000000000
23 1 17.759500000000 4.660660048033 2.570500000000
24 1 17.759500000000 2.796396028820 0.000000000000
25 1 0.000000000000 7.457056076853 2.570500000000
26 1 0.000000000000 5.592792057640 0.000000000000
27 1 1.614500000000 10.253452105673 2.570500000000
28 1 1.614500000000 8.389188086460 0.000000000000
29 1 3.229000000000 7.457056076853 2.570500000000
30 1 3.229000000000 5.592792057640 0.000000000000
31 1 4.843500000000 10.253452105673 2.570500000000
32 1 4.843500000000 8.389188086460 0.000000000000
33 1 6.458000000000 7.457056076853 2.570500000000
34 1 6.458000000000 5.592792057640 0.000000000000
35 1 8.072500000000 10.253452105673 2.570500000000
36 1 8.072500000000 8.389188086460 0.000000000000
37 1 9.687000000000 7.457056076853 2.570500000000
38 1 9.687000000000 5.592792057640 0.000000000000
39 1 11.301500000000 10.253452105673 2.570500000000
40 1 11.301500000000 8.389188086460 0.000000000000
41 1 12.916000000000 7.457056076853 2.570500000000
42 1 12.916000000000 5.592792057640 0.000000000000
43 1 14.530500000000 10.253452105673 2.570500000000
44 1 14.530500000000 8.389188086460 0.000000000000
45 1 16.145000000000 7.457056076853 2.570500000000
46 1 16.145000000000 5.592792057640 0.000000000000
47 1 17.759500000000 10.253452105673 2.570500000000
48 1 17.759500000000 8.389188086460 0.000000000000
49 1 0.000000000000 13.049848134493 2.570500000000
50 1 0.000000000000 11.185584115280 0.000000000000
51 1 1.614500000000 15.846244163313 2.570500000000
52 1 1.614500000000 13.981980144100 0.000000000000
53 1 3.229000000000 13.049848134493 2.570500000000
54 1 3.229000000000 11.185584115280 0.000000000000
55 1 4.843500000000 15.846244163313 2.570500000000
56 1 4.843500000000 13.981980144100 0.000000000000
57 1 6.458000000000 13.049848134493 2.570500000000
58 1 6.458000000000 11.185584115280 0.000000000000
59 1 8.072500000000 15.846244163313 2.570500000000
60 1 8.072500000000 13.981980144100 0.000000000000
61 1 9.687000000000 13.049848134493 2.570500000000
62 1 9.687000000000 11.185584115280 0.000000000000
63 1 11.301500000000 15.846244163313 2.570500000000
64 1 11.301500000000 13.981980144100 0.000000000000
65 1 12.916000000000 13.049848134493 2.570500000000
66 1 12.916000000000 11.185584115280 0.000000000000
67 1 14.530500000000 15.846244163313 2.570500000000
68 1 14.530500000000 13.981980144100 0.000000000000
69 1 16.145000000000 13.049848134493 2.570500000000
70 1 16.145000000000 11.185584115280 0.000000000000
71 1 17.759500000000 15.846244163313 2.570500000000
72 1 17.759500000000 13.981980144100 0.000000000000
73 1 0.000000000000 18.642640192133 2.570500000000
74 1 0.000000000000 16.778376172920 0.000000000000
75 1 1.614500000000 21.439036220953 2.570500000000
76 1 1.614500000000 19.574772201740 0.000000000000
77 1 3.229000000000 18.642640192133 2.570500000000
78 1 3.229000000000 16.778376172920 0.000000000000
79 1 4.843500000000 21.439036220953 2.570500000000
80 1 4.843500000000 19.574772201740 0.000000000000
81 1 6.458000000000 18.642640192133 2.570500000000
82 1 6.458000000000 16.778376172920 0.000000000000
83 1 8.072500000000 21.439036220953 2.570500000000
84 1 8.072500000000 19.574772201740 0.000000000000
85 1 9.687000000000 18.642640192133 2.570500000000
86 1 9.687000000000 16.778376172920 0.000000000000
87 1 11.301500000000 21.439036220953 2.570500000000
88 1 11.301500000000 19.574772201740 0.000000000000
89 1 12.916000000000 18.642640192133 2.570500000000
90 1 12.916000000000 16.778376172920 0.000000000000
91 1 14.530500000000 21.439036220953 2.570500000000
92 1 14.530500000000 19.574772201740 0.000000000000
93 1 16.145000000000 18.642640192133 2.570500000000
94 1 16.145000000000 16.778376172920 0.000000000000
95 1 17.759500000000 21.439036220953 2.570500000000
96 1 17.759500000000 19.574772201740 0.000000000000
97 1 0.000000000000 24.235432249773 2.570500000000
98 1 0.000000000000 22.371168230560 0.000000000000
99 1 1.614500000000 27.031828278593 2.570500000000
100 1 1.614500000000 25.167564259380 0.000000000000
101 1 3.229000000000 24.235432249773 2.570500000000
102 1 3.229000000000 22.371168230560 0.000000000000
103 1 4.843500000000 27.031828278593 2.570500000000
104 1 4.843500000000 25.167564259380 0.000000000000
105 1 6.458000000000 24.235432249773 2.570500000000
106 1 6.458000000000 22.371168230560 0.000000000000
107 1 8.072500000000 27.031828278593 2.570500000000
108 1 8.072500000000 25.167564259380 0.000000000000
109 1 9.687000000000 24.235432249773 2.570500000000
110 1 9.687000000000 22.371168230560 0.000000000000
111 1 11.301500000000 27.031828278593 2.570500000000
112 1 11.301500000000 25.167564259380 0.000000000000
113 1 12.916000000000 24.235432249773 2.570500000000
114 1 12.916000000000 22.371168230560 0.000000000000
115 1 14.530500000000 27.031828278593 2.570500000000
116 1 14.530500000000 25.167564259380 0.000000000000
117 1 16.145000000000 24.235432249773 2.570500000000
118 1 16.145000000000 22.371168230560 0.000000000000
119 1 17.759500000000 27.031828278593 2.570500000000
120 1 17.759500000000 25.167564259380 0.000000000000
121 1 0.000000000000 29.828224307413 2.570500000000
122 1 0.000000000000 27.963960288200 0.000000000000
123 1 1.614500000000 32.624620336233 2.570500000000
124 1 1.614500000000 30.760356317019 0.000000000000
125 1 3.229000000000 29.828224307413 2.570500000000
126 1 3.229000000000 27.963960288200 0.000000000000
127 1 4.843500000000 32.624620336233 2.570500000000
128 1 4.843500000000 30.760356317019 0.000000000000
129 1 6.458000000000 29.828224307413 2.570500000000
130 1 6.458000000000 27.963960288200 0.000000000000
131 1 8.072500000000 32.624620336233 2.570500000000
132 1 8.072500000000 30.760356317019 0.000000000000
133 1 9.687000000000 29.828224307413 2.570500000000
134 1 9.687000000000 27.963960288200 0.000000000000
135 1 11.301500000000 32.624620336233 2.570500000000
136 1 11.301500000000 30.760356317019 0.000000000000
137 1 12.916000000000 29.828224307413 2.570500000000
138 1 12.916000000000 27.963960288200 0.000000000000
139 1 14.530500000000 32.624620336233 2.570500000000
140 1 14.530500000000 30.760356317019 0.000000000000
141 1 16.145000000000 29.828224307413 2.570500000000
142 1 16.145000000000 27.963960288200 0.000000000000
143 1 17.759500000000 32.624620336233 2.570500000000
144 1 17.759500000000 30.760356317019 0.000000000000
145 1 0.000000000000 1.864264019213 7.711500000000
146 1 0.000000000000 0.000000000000 5.141000000000
147 1 1.614500000000 4.660660048033 7.711500000000
148 1 1.614500000000 2.796396028820 5.141000000000
149 1 3.229000000000 1.864264019213 7.711500000000
150 1 3.229000000000 0.000000000000 5.141000000000
151 1 4.843500000000 4.660660048033 7.711500000000
152 1 4.843500000000 2.796396028820 5.141000000000
153 1 6.458000000000 1.864264019213 7.711500000000
154 1 6.458000000000 0.000000000000 5.141000000000
155 1 8.072500000000 4.660660048033 7.711500000000
156 1 8.072500000000 2.796396028820 5.141000000000
157 1 9.687000000000 1.864264019213 7.711500000000
158 1 9.687000000000 0.000000000000 5.141000000000
159 1 11.301500000000 4.660660048033 7.711500000000
160 1 11.301500000000 2.796396028820 5.141000000000
161 1 12.916000000000 1.864264019213 7.711500000000
162 1 12.916000000000 0.000000000000 5.141000000000
163 1 14.530500000000 4.660660048033 7.711500000000
164 1 14.530500000000 2.796396028820 5.141000000000
165 1 16.145000000000 1.864264019213 7.711500000000
166 1 16.145000000000 0.000000000000 5.141000000000
167 1 17.759500000000 4.660660048033 7.711500000000
168 1 17.759500000000 2.796396028820 5.141000000000
169 1 0.000000000000 7.457056076853 7.711500000000
170 1 0.000000000000 5.592792057640 5.141000000000
171 1 1.614500000000 10.253452105673 7.711500000000
172 1 1.614500000000 8.389188086460 5.141000000000
173 1 3.229000000000 7.457056076853 7.711500000000
174 1 3.229000000000 5.592792057640 5.141000000000
175 1 4.843500000000 10.253452105673 7.711500000000
176 1 4.843500000000 8.389188086460 5.141000000000
177 1 6.458000000000 7.457056076853 7.711500000000
178 1 6.458000000000 5.592792057640 5.141000000000
179 1 8.072500000000 10.253452105673 7.711500000000
180 1 8.072500000000 8.389188086460 5.141000000000
181 1 9.687000000000 7.457056076853 7.711500000000
182 1 9.687000000000 5.592792057640 5.141000000000
183 1 11.301500000000 10.253452105673 7.711500000000
184 1 11.301500000000 8.389188086460 5.141000000000
185 1 12.916000000000 7.457056076853 7.711500000000
186 1 12.916000000000 5.592792057640 5.141000000000
187 1 14.530500000000 10.253452105673 7.711500000000
188 1 14.530500000000 8.389188086460 5.141000000000
189 1 16.145000000000 7.457056076853 7.711500000000
190 1 16.145000000000 5.592792057640 5.141000000000
191 1 17.759500000000 10.253452105673 7.711500000000
192 1 17.759500000000 8.389188086460 5.141000000000
193 1 0.000000000000 13.049848134493 7.711500000000
194 1 0.000000000000 11.185584115280 5.141000000000
195 1 1.614500000000 15.846244163313 7.711500000000
196 1 1.614500000000 13.981980144100 5.141000000000
197 1 3.229000000000 13.049848134493 7.711500000000
198 1 3.229000000000 11.185584115280 5.141000000000
199 1 4.843500000000 15.846244163313 7.711500000000
200 1 4.843500000000 13.981980144100 5.141000000000
201 1 6.458000000000 13.049848134493 7.711500000000
202 1 6.458000000000 11.185584115280 5.141000000000
203 1 8.072500000000 15.846244163313 7.711500000000
204 1 8.072500000000 13.981980144100 5.141000000000
205 1 9.687000000000 13.049848134493 7.711500000000
206 1 9.687000000000 11.185584115280 5.141000000000
207 1 11.301500000000 15.846244163313 7.711500000000
208 1 11.301500000000 13.981980144100 5.141000000000
209 1 12.916000000000 13.049848134493 7.711500000000
210 1 12.916000000000 11.185584115280 5.141000000000
211 1 14.530500000000 15.846244163313 7.711500000000
212 1 14.530500000000 13.981980144100 5.141000000000
213 1 16.145000000000 13.049848134493 7.711500000000
214 1 16.145000000000 11.185584115280 5.141000000000
215 1 17.759500000000 15.846244163313 7.711500000000
216 1 17.759500000000 13.981980144100 5.141000000000
217 1 0.000000000000 18.642640192133 7.711500000000
218 1 0.000000000000 16.778376172920 5.141000000000
219 1 1.614500000000 21.439036220953 7.711500000000
220 1 1.614500000000 19.574772201740 5.141000000000
221 1 3.229000000000 18.642640192133 7.711500000000
222 1 3.229000000000 16.778376172920 5.141000000000
223 1 4.843500000000 21.439036220953 7.711500000000
224 1 4.843500000000 19.574772201740 5.141000000000
225 1 6.458000000000 18.642640192133 7.711500000000
226 1 6.458000000000 16.778376172920 5.141000000000
227 1 8.072500000000 21.439036220953 7.711500000000
228 1 8.072500000000 19.574772201740 5.141000000000
229 1 9.687000000000 18.642640192133 7.711500000000
230 1 9.687000000000 16.778376172920 5.141000000000
231 1 11.301500000000 21.439036220953 7.711500000000
232 1 11.301500000000 19.574772201740 5.141000000000
233 1 12.916000000000 18.642640192133 7.711500000000
234 1 12.916000000000 16.778376172920 5.141000000000
235 1 14.530500000000 21.439036220953 7.711500000000
236 1 14.530500000000 19.574772201740 5.141000000000
237 1 16.145000000000 18.642640192133 7.711500000000
238 1 16.145000000000 16.778376172920 5.141000000000
239 1 17.759500000000 21.439036220953 7.711500000000
240 1 17.759500000000 19.574772201740 5.141000000000
241 1 0.000000000000 24.235432249773 7.711500000000
242 1 0.000000000000 22.371168230560 5.141000000000
243 1 1.614500000000 27.031828278593 7.711500000000
244 1 1.614500000000 25.167564259380 5.141000000000
245 1 3.229000000000 24.235432249773 7.711500000000
246 1 3.229000000000 22.371168230560 5.141000000000
247 1 4.843500000000 27.031828278593 7.711500000000
248 1 4.843500000000 25.167564259380 5.141000000000
249 1 6.458000000000 24.235432249773 7.711500000000
250 1 6.458000000000 22.371168230560 5.141000000000
251 1 8.072500000000 27.031828278593 7.711500000000
252 1 8.072500000000 25.167564259380 5.141000000000
253 1 9.687000000000 24.235432249773 7.711500000000
254 1 9.687000000000 22.371168230560 5.141000000000
255 1 11.301500000000 27.031828278593 7.711500000000
256 1 11.301500000000 25.167564259380 5.141000000000
257 1 12.916000000000 24.235432249773 7.711500000000
258 1 12.916000000000 22.371168230560 5.141000000000
259 1 14.530500000000 27.031828278593 7.711500000000
260 1 14.530500000000 25.167564259380 5.141000000000
261 1 16.145000000000 24.235432249773 7.711500000000
262 1 16.145000000000 22.371168230560 5.141000000000
263 1 17.759500000000 27.031828278593 7.711500000000
264 1 17.759500000000 25.167564259380 5.141000000000
265 1 0.000000000000 29.828224307413 7.711500000000
266 1 0.000000000000 27.963960288200 5.141000000000
267 1 1.614500000000 32.624620336233 7.711500000000
268 1 1.614500000000 30.760356317019 5.141000000000
269 1 3.229000000000 29.828224307413 7.711500000000
270 1 3.229000000000 27.963960288200 5.141000000000
271 1 4.843500000000 32.624620336233 7.711500000000
272 1 4.843500000000 30.760356317019 5.141000000000
273 1 6.458000000000 29.828224307413 7.711500000000
274 1 6.458000000000 27.963960288200 5.141000000000
275 1 8.072500000000 32.624620336233 7.711500000000
276 1 8.072500000000 30.760356317019 5.141000000000
277 1 9.687000000000 29.828224307413 7.711500000000
278 1 9.687000000000 27.963960288200 5.141000000000
279 1 11.301500000000 32.624620336233 7.711500000000
280 1 11.301500000000 30.760356317019 5.141000000000
281 1 12.916000000000 29.828224307413 7.711500000000
282 1 12.916000000000 27.963960288200 5.141000000000
283 1 14.530500000000 32.624620336233 7.711500000000
284 1 14.530500000000 30.760356317019 5.141000000000
285 1 16.145000000000 29.828224307413 7.711500000000
286 1 16.145000000000 27.963960288200 5.141000000000
287 1 17.759500000000 32.624620336233 7.711500000000
288 1 17.759500000000 30.760356317019 5.141000000000
289 1 0.000000000000 1.864264019213 12.852500000000
290 1 0.000000000000 0.000000000000 10.282000000000
291 1 1.614500000000 4.660660048033 12.852500000000
292 1 1.614500000000 2.796396028820 10.282000000000
293 1 3.229000000000 1.864264019213 12.852500000000
294 1 3.229000000000 0.000000000000 10.282000000000
295 1 4.843500000000 4.660660048033 12.852500000000
296 1 4.843500000000 2.796396028820 10.282000000000
297 1 6.458000000000 1.864264019213 12.852500000000
298 1 6.458000000000 0.000000000000 10.282000000000
299 1 8.072500000000 4.660660048033 12.852500000000
300 1 8.072500000000 2.796396028820 10.282000000000
301 1 9.687000000000 1.864264019213 12.852500000000
302 1 9.687000000000 0.000000000000 10.282000000000
303 1 11.301500000000 4.660660048033 12.852500000000
304 1 11.301500000000 2.796396028820 10.282000000000
305 1 12.916000000000 1.864264019213 12.852500000000
306 1 12.916000000000 0.000000000000 10.282000000000
307 1 14.530500000000 4.660660048033 12.852500000000
308 1 14.530500000000 2.796396028820 10.282000000000
309 1 16.145000000000 1.864264019213 12.852500000000
310 1 16.145000000000 0.000000000000 10.282000000000
311 1 17.759500000000 4.660660048033 12.852500000000
312 1 17.759500000000 2.796396028820 10.282000000000
313 1 0.000000000000 7.457056076853 12.852500000000
314 1 0.000000000000 5.592792057640 10.282000000000
315 1 1.614500000000 10.253452105673 12.852500000000
316 1 1.614500000000 8.389188086460 10.282000000000
317 1 3.229000000000 7.457056076853 12.852500000000
318 1 3.229000000000 5.592792057640 10.282000000000
319 1 4.843500000000 10.253452105673 12.852500000000
320 1 4.843500000000 8.389188086460 10.282000000000
321 1 6.458000000000 7.457056076853 12.852500000000
322 1 6.458000000000 5.592792057640 10.282000000000
323 1 8.072500000000 10.253452105673 12.852500000000
324 1 8.072500000000 8.389188086460 10.282000000000
325 1 9.687000000000 7.457056076853 12.852500000000
326 1 9.687000000000 5.592792057640 10.282000000000
327 1 11.301500000000 10.253452105673 12.852500000000
328 1 11.301500000000 8.389188086460 10.282000000000
329 1 12.916000000000 7.457056076853 12.852500000000
330 1 12.916000000000 5.592792057640 10.282000000000
331 1 14.530500000000 10.253452105673 12.852500000000
332 1 14.530500000000 8.389188086460 10.282000000000
333 1 16.145000000000 7.457056076853 12.852500000000
334 1 16.145000000000 5.592792057640 10.282000000000
335 1 17.759500000000 10.253452105673 12.852500000000
336 1 17.759500000000 8.389188086460 10.282000000000
337 1 0.000000000000 13.049848134493 12.852500000000
338 1 0.000000000000 11.185584115280 10.282000000000
339 1 1.614500000000 15.846244163313 12.852500000000
340 1 1.614500000000 13.981980144100 10.282000000000
341 1 3.229000000000 13.049848134493 12.852500000000
342 1 3.229000000000 11.185584115280 10.282000000000
343 1 4.843500000000 15.846244163313 12.852500000000
344 1 4.843500000000 13.981980144100 10.282000000000
345 1 6.458000000000 13.049848134493 12.852500000000
346 1 6.458000000000 11.185584115280 10.282000000000
347 1 8.072500000000 15.846244163313 12.852500000000
348 1 8.072500000000 13.981980144100 10.282000000000
349 1 9.687000000000 13.049848134493 12.852500000000
350 1 9.687000000000 11.185584115280 10.282000000000
351 1 11.301500000000 15.846244163313 12.852500000000
352 1 11.301500000000 13.981980144100 10.282000000000
353 1 12.916000000000 13.049848134493 12.852500000000
354 1 12.916000000000 11.185584115280 10.282000000000
355 1 14.530500000000 15.846244163313 12.852500000000
356 1 14.530500000000 13.981980144100 10.282000000000
357 1 16.145000000000 13.049848134493 12.852500000000
358 1 16.145000000000 11.185584115280 10.282000000000
359 1 17.759500000000 15.846244163313 12.852500000000
360 1 17.759500000000 13.981980144100 10.282000000000
361 1 0.000000000000 18.642640192133 12.852500000000
362 1 0.000000000000 16.778376172920 10.282000000000
363 1 1.614500000000 21.439036220953 12.852500000000
364 1 1.614500000000 19.574772201740 10.282000000000
365 1 3.229000000000 18.642640192133 12.852500000000
366 1 3.229000000000 16.778376172920 10.282000000000
367 1 4.843500000000 21.439036220953 12.852500000000
368 1 4.843500000000 19.574772201740 10.282000000000
369 1 6.458000000000 18.642640192133 12.852500000000
370 1 6.458000000000 16.778376172920 10.282000000000
371 1 8.072500000000 21.439036220953 12.852500000000
372 1 8.072500000000 19.574772201740 10.282000000000
373 1 9.687000000000 18.642640192133 12.852500000000
374 1 9.687000000000 16.778376172920 10.282000000000
375 1 11.301500000000 21.439036220953 12.852500000000
376 1 11.301500000000 19.574772201740 10.282000000000
377 1 12.916000000000 18.642640192133 12.852500000000
378 1 12.916000000000 16.778376172920 10.282000000000
379 1 14.530500000000 21.439036220953 12.852500000000
380 1 14.530500000000 19.574772201740 10.282000000000
381 1 16.145000000000 18.642640192133 12.852500000000
382 1 16.145000000000 16.778376172920 10.282000000000
383 1 17.759500000000 21.439036220953 12.852500000000
384 1 17.759500000000 19.574772201740 10.282000000000
385 1 0.000000000000 24.235432249773 12.852500000000
386 1 0.000000000000 22.371168230560 10.282000000000
387 1 1.614500000000 27.031828278593 12.852500000000
388 1 1.614500000000 25.167564259380 10.282000000000
389 1 3.229000000000 24.235432249773 12.852500000000
390 1 3.229000000000 22.371168230560 10.282000000000
391 1 4.843500000000 27.031828278593 12.852500000000
392 1 4.843500000000 25.167564259380 10.282000000000
393 1 6.458000000000 24.235432249773 12.852500000000
394 1 6.458000000000 22.371168230560 10.282000000000
395 1 8.072500000000 27.031828278593 12.852500000000
396 1 8.072500000000 25.167564259380 10.282000000000
397 1 9.687000000000 24.235432249773 12.852500000000
398 1 9.687000000000 22.371168230560 10.282000000000
399 1 11.301500000000 27.031828278593 12.852500000000
400 1 11.301500000000 25.167564259380 10.282000000000
401 1 12.916000000000 24.235432249773 12.852500000000
402 1 12.916000000000 22.371168230560 10.282000000000
403 1 14.530500000000 27.031828278593 12.852500000000
404 1 14.530500000000 25.167564259380 10.282000000000
405 1 16.145000000000 24.235432249773 12.852500000000
406 1 16.145000000000 22.371168230560 10.282000000000
407 1 17.759500000000 27.031828278593 12.852500000000
408 1 17.759500000000 25.167564259380 10.282000000000
409 1 0.000000000000 29.828224307413 12.852500000000
410 1 0.000000000000 27.963960288200 10.282000000000
411 1 1.614500000000 32.624620336233 12.852500000000
412 1 1.614500000000 30.760356317019 10.282000000000
413 1 3.229000000000 29.828224307413 12.852500000000
414 1 3.229000000000 27.963960288200 10.282000000000
415 1 4.843500000000 32.624620336233 12.852500000000
416 1 4.843500000000 30.760356317019 10.282000000000
417 1 6.458000000000 29.828224307413 12.852500000000
418 1 6.458000000000 27.963960288200 10.282000000000
419 1 8.072500000000 32.624620336233 12.852500000000
420 1 8.072500000000 30.760356317019 10.282000000000
421 1 9.687000000000 29.828224307413 12.852500000000
422 1 9.687000000000 27.963960288200 10.282000000000
423 1 11.301500000000 32.624620336233 12.852500000000
424 1 11.301500000000 30.760356317019 10.282000000000
425 1 12.916000000000 29.828224307413 12.852500000000
426 1 12.916000000000 27.963960288200 10.282000000000
427 1 14.530500000000 32.624620336233 12.852500000000
428 1 14.530500000000 30.760356317019 10.282000000000
429 1 16.145000000000 29.828224307413 12.852500000000
430 1 16.145000000000 27.963960288200 10.282000000000
431 1 17.759500000000 32.624620336233 12.852500000000
432 1 17.759500000000 30.760356317019 10.282000000000
433 1 0.000000000000 1.864264019213 17.993500000000
434 1 0.000000000000 0.000000000000 15.423000000000
435 1 1.614500000000 4.660660048033 17.993500000000
436 1 1.614500000000 2.796396028820 15.423000000000
437 1 3.229000000000 1.864264019213 17.993500000000
438 1 3.229000000000 0.000000000000 15.423000000000
439 1 4.843500000000 4.660660048033 17.993500000000
440 1 4.843500000000 2.796396028820 15.423000000000
441 1 6.458000000000 1.864264019213 17.993500000000
442 1 6.458000000000 0.000000000000 15.423000000000
443 1 8.072500000000 4.660660048033 17.993500000000
444 1 8.072500000000 2.796396028820 15.423000000000
445 1 9.687000000000 1.864264019213 17.993500000000
446 1 9.687000000000 0.000000000000 15.423000000000
447 1 11.301500000000 4.660660048033 17.993500000000
448 1 11.301500000000 2.796396028820 15.423000000000
449 1 12.916000000000 1.864264019213 17.993500000000
450 1 12.916000000000 0.000000000000 15.423000000000
451 1 14.530500000000 4.660660048033 17.993500000000
452 1 14.530500000000 2.796396028820 15.423000000000
453 1 16.145000000000 1.864264019213 17.993500000000
454 1 16.145000000000 0.000000000000 15.423000000000
455 1 17.759500000000 4.660660048033 17.993500000000
456 1 17.759500000000 2.796396028820 15.423000000000
457 1 0.000000000000 7.457056076853 17.993500000000
458 1 0.000000000000 5.592792057640 15.423000000000
459 1 1.614500000000 10.253452105673 17.993500000000
460 1 1.614500000000 8.389188086460 15.423000000000
461 1 3.229000000000 7.457056076853 17.993500000000
462 1 3.229000000000 5.592792057640 15.423000000000
463 1 4.843500000000 10.253452105673 17.993500000000
464 1 4.843500000000 8.389188086460 15.423000000000
465 1 6.458000000000 7.457056076853 17.993500000000
466 1 6.458000000000 5.592792057640 15.423000000000
467 1 8.072500000000 10.253452105673 17.993500000000
468 1 8.072500000000 8.389188086460 15.423000000000
469 1 9.687000000000 7.457056076853 17.993500000000
470 1 9.687000000000 5.592792057640 15.423000000000
471 1 11.301500000000 10.253452105673 17.993500000000
472 1 11.301500000000 8.389188086460 15.423000000000
473 1 12.916000000000 7.457056076853 17.993500000000
474 1 12.916000000000 5.592792057640 15.423000000000
475 1 14.530500000000 10.253452105673 17.993500000000
476 1 14.530500000000 8.389188086460 15.423000000000
477 1 16.145000000000 7.457056076853 17.993500000000
478 1 16.145000000000 5.592792057640 15.423000000000
479 1 17.759500000000 10.253452105673 17.993500000000
480 1 17.759500000000 8.389188086460 15.423000000000
481 1 0.000000000000 13.049848134493 17.993500000000
482 1 0.000000000000 11.185584115280 15.423000000000
483 1 1.614500000000 15.846244163313 17.993500000000
484 1 1.614500000000 13.981980144100 15.423000000000
485 1 3.229000000000 13.049848134493 17.993500000000
486 1 3.229000000000 11.185584115280 15.423000000000
487 1 4.843500000000 15.846244163313 17.993500000000
488 1 4.843500000000 13.981980144100 15.423000000000
489 1 6.458000000000 13.049848134493 17.993500000000
490 1 6.458000000000 11.185584115280 15.423000000000
491 1 8.072500000000 15.846244163313 17.993500000000
492 1 8.072500000000 13.981980144100 15.423000000000
493 1 9.687000000000 13.049848134493 17.993500000000
494 1 9.687000000000 11.185584115280 15.423000000000
495 1 11.301500000000 15.846244163313 17.993500000000
496 1 11.301500000000 13.981980144100 15.423000000000
497 1 12.916000000000 13.049848134493 17.993500000000
498 1 12.916000000000 11.185584115280 15.423000000000
499 1 14.530500000000 15.846244163313 17.993500000000
500 1 14.530500000000 13.981980144100 15.423000000000
501 1 16.145000000000 13.049848134493 17.993500000000
502 1 16.145000000000 11.185584115280 15.423000000000
503 1 17.759500000000 15.846244163313 17.993500000000
504 1 17.759500000000 13.981980144100 15.423000000000
505 1 0.000000000000 18.642640192133 17.993500000000
506 1 0.000000000000 16.778376172920 15.423000000000
507 1 1.614500000000 21.439036220953 17.993500000000
508 1 1.614500000000 19.574772201740 15.423000000000
509 1 3.229000000000 18.642640192133 17.993500000000
510 1 3.229000000000 16.778376172920 15.423000000000
511 1 4.843500000000 21.439036220953 17.993500000000
512 1 4.843500000000 19.574772201740 15.423000000000
513 1 6.458000000000 18.642640192133 17.993500000000
514 1 6.458000000000 16.778376172920 15.423000000000
515 1 8.072500000000 21.439036220953 17.993500000000
516 1 8.072500000000 19.574772201740 15.423000000000
517 1 9.687000000000 18.642640192133 17.993500000000
518 1 9.687000000000 16.778376172920 15.423000000000
519 1 11.301500000000 21.439036220953 17.993500000000
520 1 11.301500000000 19.574772201740 15.423000000000
521 1 12.916000000000 18.642640192133 17.993500000000
522 1 12.916000000000 16.778376172920 15.423000000000
523 1 14.530500000000 21.439036220953 17.993500000000
524 1 14.530500000000 19.574772201740 15.423000000000
525 1 16.145000000000 18.642640192133 17.993500000000
526 1 16.145000000000 16.778376172920 15.423000000000
527 1 17.759500000000 21.439036220953 17.993500000000
528 1 17.759500000000 19.574772201740 15.423000000000
529 1 0.000000000000 24.235432249773 17.993500000000
530 1 0.000000000000 22.371168230560 15.423000000000
531 1 1.614500000000 27.031828278593 17.993500000000
532 1 1.614500000000 25.167564259380 15.423000000000
533 1 3.229000000000 24.235432249773 17.993500000000
534 1 3.229000000000 22.371168230560 15.423000000000
535 1 4.843500000000 27.031828278593 17.993500000000
536 1 4.843500000000 25.167564259380 15.423000000000
537 1 6.458000000000 24.235432249773 17.993500000000
538 1 6.458000000000 22.371168230560 15.423000000000
539 1 8.072500000000 27.031828278593 17.993500000000
540 1 8.072500000000 25.167564259380 15.423000000000
541 1 9.687000000000 24.235432249773 17.993500000000
542 1 9.687000000000 22.371168230560 15.423000000000
543 1 11.301500000000 27.031828278593 17.993500000000
544 1 11.301500000000 25.167564259380 15.423000000000
545 1 12.916000000000 24.235432249773 17.993500000000
546 1 12.916000000000 22.371168230560 15.423000000000
547 1 14.530500000000 27.031828278593 17.993500000000
548 1 14.530500000000 25.167564259380 15.423000000000
549 1 16.145000000000 24.235432249773 17.993500000000
550 1 16.145000000000 22.371168230560 15.423000000000
551 1 17.759500000000 27.031828278593 17.993500000000
552 1 17.759500000000 25.167564259380 15.423000000000
553 1 0.000000000000 29.828224307413 17.993500000000
554 1 0.000000000000 27.963960288200 15.423000000000
555 1 1.614500000000 32.624620336233 17.993500000000
556 1 1.614500000000 30.760356317019 15.423000000000
557 1 3.229000000000 29.828224307413 17.993500000000
558 1 3.229000000000 27.963960288200 15.423000000000
559 1 4.843500000000 32.624620336233 17.993500000000
560 1 4.843500000000 30.760356317019 15.423000000000
561 1 6.458000000000 29.828224307413 17.993500000000
562 1 6.458000000000 27.963960288200 15.423000000000
563 1 8.072500000000 32.624620336233 17.993500000000
564 1 8.072500000000 30.760356317019 15.423000000000
565 1 9.687000000000 29.828224307413 17.993500000000
566 1 9.687000000000 27.963960288200 15.423000000000
567 1 11.301500000000 32.624620336233 17.993500000000
568 1 11.301500000000 30.760356317019 15.423000000000
569 1 12.916000000000 29.828224307413 17.993500000000
570 1 12.916000000000 27.963960288200 15.423000000000
571 1 14.530500000000 32.624620336233 17.993500000000
572 1 14.530500000000 30.760356317019 15.423000000000
573 1 16.145000000000 29.828224307413 17.993500000000
574 1 16.145000000000 27.963960288200 15.423000000000
575 1 17.759500000000 32.624620336233 17.993500000000
576 1 17.759500000000 30.760356317019 15.423000000000
577 1 0.000000000000 1.864264019213 23.134500000000
578 1 0.000000000000 0.000000000000 20.564000000000
579 1 1.614500000000 4.660660048033 23.134500000000
580 1 1.614500000000 2.796396028820 20.564000000000
581 1 3.229000000000 1.864264019213 23.134500000000
582 1 3.229000000000 0.000000000000 20.564000000000
583 1 4.843500000000 4.660660048033 23.134500000000
584 1 4.843500000000 2.796396028820 20.564000000000
585 1 6.458000000000 1.864264019213 23.134500000000
586 1 6.458000000000 0.000000000000 20.564000000000
587 1 8.072500000000 4.660660048033 23.134500000000
588 1 8.072500000000 2.796396028820 20.564000000000
589 1 9.687000000000 1.864264019213 23.134500000000
590 1 9.687000000000 0.000000000000 20.564000000000
591 1 11.301500000000 4.660660048033 23.134500000000
592 1 11.301500000000 2.796396028820 20.564000000000
593 1 12.916000000000 1.864264019213 23.134500000000
594 1 12.916000000000 0.000000000000 20.564000000000
595 1 14.530500000000 4.660660048033 23.134500000000
596 1 14.530500000000 2.796396028820 20.564000000000
597 1 16.145000000000 1.864264019213 23.134500000000
598 1 16.145000000000 0.000000000000 20.564000000000
599 1 17.759500000000 4.660660048033 23.134500000000
600 1 17.759500000000 2.796396028820 20.564000000000
601 1 0.000000000000 7.457056076853 23.134500000000
602 1 0.000000000000 5.592792057640 20.564000000000
603 1 1.614500000000 10.253452105673 23.134500000000
604 1 1.614500000000 8.389188086460 20.564000000000
605 1 3.229000000000 7.457056076853 23.134500000000
606 1 3.229000000000 5.592792057640 20.564000000000
607 1 4.843500000000 10.253452105673 23.134500000000
608 1 4.843500000000 8.389188086460 20.564000000000
609 1 6.458000000000 7.457056076853 23.134500000000
610 1 6.458000000000 5.592792057640 20.564000000000
611 1 8.072500000000 10.253452105673 23.134500000000
612 1 8.072500000000 8.389188086460 20.564000000000
613 1 9.687000000000 7.457056076853 23.134500000000
614 1 9.687000000000 5.592792057640 20.564000000000
615 1 11.301500000000 10.253452105673 23.134500000000
616 1 11.301500000000 8.389188086460 20.564000000000
617 1 12.916000000000 7.457056076853 23.134500000000
618 1 12.916000000000 5.592792057640 20.564000000000
619 1 14.530500000000 10.253452105673 23.134500000000
620 1 14.530500000000 8.389188086460 20.564000000000
621 1 16.145000000000 7.457056076853 23.134500000000
622 1 16.145000000000 5.592792057640 20.564000000000
623 1 17.759500000000 10.253452105673 23.134500000000
624 1 17.759500000000 8.389188086460 20.564000000000
625 1 0.000000000000 13.049848134493 23.134500000000
626 1 0.000000000000 11.185584115280 20.564000000000
627 1 1.614500000000 15.846244163313 23.134500000000
628 1 1.614500000000 13.981980144100 20.564000000000
629 1 3.229000000000 13.049848134493 23.134500000000
630 1 3.229000000000 11.185584115280 20.564000000000
631 1 4.843500000000 15.846244163313 23.134500000000
632 1 4.843500000000 13.981980144100 20.564000000000
633 1 6.458000000000 13.049848134493 23.134500000000
634 1 6.458000000000 11.185584115280 20.564000000000
635 1 8.072500000000 15.846244163313 23.134500000000
636 1 8.072500000000 13.981980144100 20.564000000000
637 1 9.687000000000 13.049848134493 23.134500000000
638 1 9.687000000000 11.185584115280 20.564000000000
639 1 11.301500000000 15.846244163313 23.134500000000
640 1 11.301500000000 13.981980144100 20.564000000000
641 1 12.916000000000 13.049848134493 23.134500000000
642 1 12.916000000000 11.185584115280 20.564000000000
643 1 14.530500000000 15.846244163313 23.134500000000
644 1 14.530500000000 13.981980144100 20.564000000000
645 1 16.145000000000 13.049848134493 23.134500000000
646 1 16.145000000000 11.185584115280 20.564000000000
647 1 17.759500000000 15.846244163313 23.134500000000
648 1 17.759500000000 13.981980144100 20.564000000000
649 1 0.000000000000 18.642640192133 23.134500000000
650 1 0.000000000000 16.778376172920 20.564000000000
651 1 1.614500000000 21.439036220953 23.134500000000
652 1 1.614500000000 19.574772201740 20.564000000000
653 1 3.229000000000 18.642640192133 23.134500000000
654 1 3.229000000000 16.778376172920 20.564000000000
655 1 4.843500000000 21.439036220953 23.134500000000
656 1 4.843500000000 19.574772201740 20.564000000000
657 1 6.458000000000 18.642640192133 23.134500000000
658 1 6.458000000000 16.778376172920 20.564000000000
659 1 8.072500000000 21.439036220953 23.134500000000
660 1 8.072500000000 19.574772201740 20.564000000000
661 1 9.687000000000 18.642640192133 23.134500000000
662 1 9.687000000000 16.778376172920 20.564000000000
663 1 11.301500000000 21.439036220953 23.134500000000
664 1 11.301500000000 19.574772201740 20.564000000000
665 1 12.916000000000 18.642640192133 23.134500000000
666 1 12.916000000000 16.778376172920 20.564000000000
667 1 14.530500000000 21.439036220953 23.134500000000
668 1 14.530500000000 19.574772201740 20.564000000000
669 1 16.145000000000 18.642640192133 23.134500000000
670 1 16.145000000000 16.778376172920 20.564000000000
671 1 17.759500000000 21.439036220953 23.134500000000
672 1 17.759500000000 19.574772201740 20.564000000000
673 1 0.000000000000 24.235432249773 23.134500000000
674 1 0.000000000000 22.371168230560 20.564000000000
675 1 1.614500000000 27.031828278593 23.134500000000
676 1 1.614500000000 25.167564259380 20.564000000000
677 1 3.229000000000 24.235432249773 23.134500000000
678 1 3.229000000000 22.371168230560 20.564000000000
679 1 4.843500000000 27.031828278593 23.134500000000
680 1 4.843500000000 25.167564259380 20.564000000000
681 1 6.458000000000 24.235432249773 23.134500000000
682 1 6.458000000000 22.371168230560 20.564000000000
683 1 8.072500000000 27.031828278593 23.134500000000
684 1 8.072500000000 25.167564259380 20.564000000000
685 1 9.687000000000 24.235432249773 23.134500000000
686 1 9.687000000000 22.371168230560 20.564000000000
687 1 11.301500000000 27.031828278593 23.134500000000
688 1 11.301500000000 25.167564259380 20.564000000000
689 1 12.916000000000 24.235432249773 23.134500000000
690 1 12.916000000000 22.371168230560 20.564000000000
691 1 14.530500000000 27.031828278593 23.134500000000
692 1 14.530500000000 25.167564259380 20.564000000000
693 1 16.145000000000 24.235432249773 23.134500000000
694 1 16.145000000000 22.371168230560 20.564000000000
695 1 17.759500000000 27.031828278593 23.134500000000
696 1 17.759500000000 25.167564259380 20.564000000000
697 1 0.000000000000 29.828224307413 23.134500000000
698 1 0.000000000000 27.963960288200 20.564000000000
699 1 1.614500000000 32.624620336233 23.134500000000
700 1 1.614500000000 30.760356317019 20.564000000000
701 1 3.229000000000 29.828224307413 23.134500000000
702 1 3.229000000000 27.963960288200 20.564000000000
703 1 4.843500000000 32.624620336233 23.134500000000
704 1 4.843500000000 30.760356317019 20.564000000000
705 1 6.458000000000 29.828224307413 23.134500000000
706 1 6.458000000000 27.963960288200 20.564000000000
707 1 8.072500000000 32.624620336233 23.134500000000
708 1 8.072500000000 30.760356317019 20.564000000000
709 1 9.687000000000 29.828224307413 23.134500000000
710 1 9.687000000000 27.963960288200 20.564000000000
711 1 11.301500000000 32.624620336233 23.134500000000
712 1 11.301500000000 30.760356317019 20.564000000000
713 1 12.916000000000 29.828224307413 23.134500000000
714 1 12.916000000000 27.963960288200 20.564000000000
715 1 14.530500000000 32.624620336233 23.134500000000
716 1 14.530500000000 30.760356317019 20.564000000000
717 1 16.145000000000 29.828224307413 23.134500000000
718 1 16.145000000000 27.963960288200 20.564000000000
719 1 17.759500000000 32.624620336233 23.134500000000
720 1 17.759500000000 30.760356317019 20.564000000000
721 1 0.000000000000 1.864264019213 28.275500000000
722 1 0.000000000000 0.000000000000 25.705000000000
723 1 1.614500000000 4.660660048033 28.275500000000
724 1 1.614500000000 2.796396028820 25.705000000000
725 1 3.229000000000 1.864264019213 28.275500000000
726 1 3.229000000000 0.000000000000 25.705000000000
727 1 4.843500000000 4.660660048033 28.275500000000
728 1 4.843500000000 2.796396028820 25.705000000000
729 1 6.458000000000 1.864264019213 28.275500000000
730 1 6.458000000000 0.000000000000 25.705000000000
731 1 8.072500000000 4.660660048033 28.275500000000
732 1 8.072500000000 2.796396028820 25.705000000000
733 1 9.687000000000 1.864264019213 28.275500000000
734 1 9.687000000000 0.000000000000 25.705000000000
735 1 11.301500000000 4.660660048033 28.275500000000
736 1 11.301500000000 2.796396028820 25.705000000000
737 1 12.916000000000 1.864264019213 28.275500000000
738 1 12.916000000000 0.000000000000 25.705000000000
739 1 14.530500000000 4.660660048033 28.275500000000
740 1 14.530500000000 2.796396028820 25.705000000000
741 1 16.145000000000 1.864264019213 28.275500000000
742 1 16.145000000000 0.000000000000 25.705000000000
743 1 17.759500000000 4.660660048033 28.275500000000
744 1 17.759500000000 2.796396028820 25.705000000000
745 1 0.000000000000 7.457056076853 28.275500000000
746 1 0.000000000000 5.592792057640 25.705000000000
747 1 1.614500000000 10.253452105673 28.275500000000
748 1 1.614500000000 8.389188086460 25.705000000000
749 1 3.229000000000 7.457056076853 28.275500000000
750 1 3.229000000000 5.592792057640 25.705000000000
751 1 4.843500000000 10.253452105673 28.275500000000
752 1 4.843500000000 8.389188086460 25.705000000000
753 1 6.458000000000 7.457056076853 28.275500000000
754 1 6.458000000000 5.592792057640 25.705000000000
755 1 8.072500000000 10.253452105673 28.275500000000
756 1 8.072500000000 8.389188086460 25.705000000000
757 1 9.687000000000 7.457056076853 28.275500000000
758 1 9.687000000000 5.592792057640 25.705000000000
759 1 11.301500000000 10.253452105673 28.275500000000
760 1 11.301500000000 8.389188086460 25.705000000000
761 1 12.916000000000 7.457056076853 28.275500000000
762 1 12.916000000000 5.592792057640 25.705000000000
763 1 14.530500000000 10.253452105673 28.275500000000
764 1 14.530500000000 8.389188086460 25.705000000000
765 1 16.145000000000 7.457056076853 28.275500000000
766 1 16.145000000000 5.592792057640 25.705000000000
767 1 17.759500000000 10.253452105673 28.275500000000
768 1 17.759500000000 8.389188086460 25.705000000000
769 1 0.000000000000 13.049848134493 28.275500000000
770 1 0.000000000000 11.185584115280 25.705000000000
771 1 1.614500000000 15.846244163313 28.275500000000
772 1 1.614500000000 13.981980144100 25.705000000000
773 1 3.229000000000 13.049848134493 28.275500000000
774 1 3.229000000000 11.185584115280 25.705000000000
775 1 4.843500000000 15.846244163313 28.275500000000
776 1 4.843500000000 13.981980144100 25.705000000000
777 1 6.458000000000 13.049848134493 28.275500000000
778 1 6.458000000000 11.185584115280 25.705000000000
779 1 8.072500000000 15.846244163313 28.275500000000
780 1 8.072500000000 13.981980144100 25.705000000000
781 1 9.687000000000 13.049848134493 28.275500000000
782 1 9.687000000000 11.185584115280 25.705000000000
783 1 11.301500000000 15.846244163313 28.275500000000
784 1 11.301500000000 13.981980144100 25.705000000000
785 1 12.916000000000 13.049848134493 28.275500000000
786 1 12.916000000000 11.185584115280 25.705000000000
787 1 14.530500000000 15.846244163313 28.275500000000
788 1 14.530500000000 13.981980144100 25.705000000000
789 1 16.145000000000 13.049848134493 28.275500000000
790 1 16.145000000000 11.185584115280 25.705000000000
791 1 17.759500000000 15.846244163313 28.275500000000
792 1 17.759500000000 13.981980144100 25.705000000000
793 1 0.000000000000 18.642640192133 28.275500000000
794 1 0.000000000000 16.778376172920 25.705000000000
795 1 1.614500000000 21.439036220953 28.275500000000
796 1 1.614500000000 19.574772201740 25.705000000000
797 1 3.229000000000 18.642640192133 28.275500000000
798 1 3.229000000000 16.778376172920 25.705000000000
799 1 4.843500000000 21.439036220953 28.275500000000
800 1 4.843500000000 19.574772201740 25.705000000000
801 1 6.458000000000 18.642640192133 28.275500000000
802 1 6.458000000000 16.778376172920 25.705000000000
803 1 8.072500000000 21.439036220953 28.275500000000
804 1 8.072500000000 19.574772201740 25.705000000000
805 1 9.687000000000 18.642640192133 28.275500000000
806 1 9.687000000000 16.778376172920 25.705000000000
807 1 11.301500000000 21.439036220953 28.275500000000
808 1 11.301500000000 19.574772201740 25.705000000000
809 1 12.916000000000 18.642640192133 28.275500000000
810 1 12.916000000000 16.778376172920 25.705000000000
811 1 14.530500000000 21.439036220953 28.275500000000
812 1 14.530500000000 19.574772201740 25.705000000000
813 1 16.145000000000 18.642640192133 28.275500000000
814 1 16.145000000000 16.778376172920 25.705000000000
815 1 17.759500000000 21.439036220953 28.275500000000
816 1 17.759500000000 19.574772201740 25.705000000000
817 1 0.000000000000 24.235432249773 28.275500000000
818 1 0.000000000000 22.371168230560 25.705000000000
819 1 1.614500000000 27.031828278593 28.275500000000
820 1 1.614500000000 25.167564259380 25.705000000000
821 1 3.229000000000 24.235432249773 28.275500000000
822 1 3.229000000000 22.371168230560 25.705000000000
823 1 4.843500000000 27.031828278593 28.275500000000
824 1 4.843500000000 25.167564259380 25.705000000000
825 1 6.458000000000 24.235432249773 28.275500000000
826 1 6.458000000000 22.371168230560 25.705000000000
827 1 8.072500000000 27.031828278593 28.275500000000
828 1 8.072500000000 25.167564259380 25.705000000000
829 1 9.687000000000 24.235432249773 28.275500000000
830 1 9.687000000000 22.371168230560 25.705000000000
831 1 11.301500000000 27.031828278593 28.275500000000
832 1 11.301500000000 25.167564259380 25.705000000000
833 1 12.916000000000 24.235432249773 28.275500000000
834 1 12.916000000000 22.371168230560 25.705000000000
835 1 14.530500000000 27.031828278593 28.275500000000
836 1 14.530500000000 25.167564259380 25.705000000000
837 1 16.145000000000 24.235432249773 28.275500000000
838 1 16.145000000000 22.371168230560 25.705000000000
839 1 17.759500000000 27.031828278593 28.275500000000
840 1 17.759500000000 25.167564259380 25.705000000000
841 1 0.000000000000 29.828224307413 28.275500000000
842 1 0.000000000000 27.963960288200 25.705000000000
843 1 1.614500000000 32.624620336233 28.275500000000
844 1 1.614500000000 30.760356317019 25.705000000000
845 1 3.229000000000 29.828224307413 28.275500000000
846 1 3.229000000000 27.963960288200 25.705000000000
847 1 4.843500000000 32.624620336233 28.275500000000
848 1 4.843500000000 30.760356317019 25.705000000000
849 1 6.458000000000 29.828224307413 28.275500000000
850 1 6.458000000000 27.963960288200 25.705000000000
851 1 8.072500000000 32.624620336233 28.275500000000
852 1 8.072500000000 30.760356317019 25.705000000000
853 1 9.687000000000 29.828224307413 28.275500000000
854 1 9.687000000000 27.963960288200 25.705000000000
855 1 11.301500000000 32.624620336233 28.275500000000
856 1 11.301500000000 30.760356317019 25.705000000000
857 1 12.916000000000 29.828224307413 28.275500000000
858 1 12.916000000000 27.963960288200 25.705000000000
859 1 14.530500000000 32.624620336233 28.275500000000
860 1 14.530500000000 30.760356317019 25.705000000000
861 1 16.145000000000 29.828224307413 28.275500000000
862 1 16.145000000000 27.963960288200 25.705000000000
863 1 17.759500000000 32.624620336233 28.275500000000
864 1 17.759500000000 30.760356317019 25.705000000000

55
examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lda_scalings.dat Normal file
View File

@ -0,0 +1,55 @@
0.65552758 -0.08218108 -0.23122826
2.03065849 0.46494117 0.87297750
-15.80180341 1.50484584 0.31669351
0.06060238 -0.47589059 -0.41017499
4.23928030 -2.27982958 4.87969884
-1.09746642 -1.28258171 2.03459312
5.48480653 -0.66345012 3.18471732
-1.57479966 0.17478998 -0.17156696
3.85779786 1.59890578 1.78936017
-1.14469715 -2.15823271 2.14353632
5.97160056 0.11573423 0.97653410
4.44645807 -0.15365582 -0.08773622
3.09452721 0.32439223 1.19779688
-1.22585061 -0.32185613 0.03949731
0.44816997 -1.11182687 0.26222208
0.19532128 0.30397832 -0.57154050
5.52432571 -0.76685448 0.32647935
6.37957282 -0.96148815 1.53439397
2.73798648 -0.69516327 1.73607004
0.94755899 0.41154702 -0.14095753
1.50733544 1.22254481 0.26284605
0.98313431 -1.24195379 0.59009611
-0.76518592 0.11605047 -0.00304658
-0.68335076 0.48935564 -0.53834507
1.86534260 -0.49032664 -0.06298849
1.52931829 0.64853878 -0.56286214
2.64217062 -1.37348638 0.22526281
0.18023516 0.03439864 0.77624538
2.02366558 0.35432524 0.76748492
0.80982907 0.31806067 0.08774175
1.57388194 -1.07822533 0.15886237
0.41345498 0.38916338 -0.29917607
-0.24819893 0.13763422 0.45471609
-2.27933523 -0.01771636 -0.20567577
1.52275665 0.35306670 0.21266257
0.28547991 1.05230832 1.16641438
0.97147437 -0.63973458 -0.37994470
0.48124764 0.03483500 -0.01982056
0.74502588 0.14367872 -0.24443596
0.48813660 0.15632903 -0.88469078
0.04886450 0.00882595 -0.47920447
0.03103900 -0.15091487 -0.41193682
-0.10106190 0.14911569 0.10727243
-0.15552036 -0.49286545 -0.04644942
0.27304084 0.35638954 1.13331445
0.57788886 -0.50269555 0.09110942
0.36780762 -0.08710371 -0.28478716
1.01678932 -0.42099561 -0.07317253
0.06561086 -0.27253002 -0.05366136
0.22266923 0.19999531 -0.30017173
-0.18666193 0.02576273 0.27752106
-0.76718071 0.61299522 0.58296511
0.60978530 0.04962900 -0.32796430
-0.11572649 0.03034386 -0.83005753
0.12675714 0.00004617 -0.37078106

1
examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_bias.dat Normal file
View File

@ -0,0 +1 @@
-6.32012657 5.62127377 1.19871662 -0.49986382

4
examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_decision.dat Normal file
View File

@ -0,0 +1,4 @@
-0.42810669 1.25467216 0.93144383
0.09624929 -0.80420088 0.48996738
-0.09865949 0.39991755 -0.69233982
0.43051689 -0.85038883 -0.72907140

20
examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mahalanobis_file.dat Normal file
View File

@ -0,0 +1,20 @@
5.0540
-23.8329 4.6638 3.9805
1.1377 0.1077 -0.0171
0.1077 0.8846 -0.2577
-0.0171 -0.2577 0.6783
5.2340
-21.2853 -6.1583 1.7948
1.7124 0.0341 0.1966
0.0341 0.6453 0.2880
0.1966 0.2880 1.8991
5.0360
-23.1593 1.3059 -5.7549
0.7496 -0.0806 -0.1101
-0.0806 1.1178 0.1667
-0.1101 0.1667 0.6711
7.9940
68.1971 0.1604 -0.0067
0.9663 -0.1846 0.6622
-0.1846 8.2371 0.9841
0.6622 0.9841 5.9601

55
examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mean_descriptor.dat Normal file
View File

@ -0,0 +1,55 @@
137.71497059
0.36342014
-2.78949838
1.75623090
-4.86893969
-2.31918628
-3.01873942
59.70217846
-8.31239311
-1.05113276
-4.08948813
11.70560234
17.48710737
42.43158755
-6.27727395
-1.46675636
-3.40739849
1.58674150
13.02515977
5.67885926
6.45692906
4.69273492
21.59764216
-7.68805780
-4.37357550
-5.79764719
0.53149261
-0.00723980
-2.47811316
-0.34939237
-4.59425510
-4.44056296
107.64051985
-9.32851480
-6.62214151
-5.69590145
22.80361437
9.47641390
2.25214024
-0.19403065
3.05386205
12.91756406
135.15381317
-9.93292065
-3.73311129
10.67039500
9.60945072
-0.03566872
21.97944941
6.70251772
74.60284853
-5.99090678
0.21877973
-1.19909174
1.37424965

57
examples/PACKAGES/sna_nnn_slcsa/in.slcsa Normal file
View File

@ -0,0 +1,57 @@
variable trequis equal 750.0
variable prequis_low equal 0.0
variable prequis_high equal 25.0e4
variable equilSteps equal 200
variable runSteps equal 2000
variable freqdump equal 200
variable pstime equal step*dt
variable sxx equal 1.e-4*pxx
variable syy equal 1.e-4*pyy
variable szz equal 1.e-4*pzz
variable sxy equal 1.e-4*pxy
variable sxz equal 1.e-4*pxz
variable syz equal 1.e-4*pyz
variable TK equal temp
variable PE equal pe
variable KE equal ke
variable V equal vol
dimension 3
boundary p p p
units metal
atom_style atomic
read_data data.zr_cell
replicate 1 5 5
change_box all triclinic
pair_style hybrid/overlay zero 9.0 eam/fs
pair_coeff * * zero
pair_coeff * * eam/fs Zr_mm.eam.fs Zr
timestep 0.002
thermo 50
thermo_style custom step pe ke temp vol pxx pyy pzz pxy pyz pxz
# fix extra all print 50 "${pstime} ${TK} ${PE} ${KE} ${V} ${sxx} ${syy} ${szz} ${sxy} ${sxz} ${syz}" file thermo_global_npt_low_temperature_Zr_hcp.dat
velocity all create ${trequis} 42345 dist gaussian
# 1st step : compute the bispectrum on 24 nearest neighbors
compute bnnn all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1 delta 0.25
# 2nd step : perform dimension reduction + logistic regression
compute slcsa all slcsa/atom 8 4 dir.slcsa/mean_descriptor.dat dir.slcsa/lda_scalings.dat dir.slcsa/lr_decision.dat dir.slcsa/lr_bias.dat dir.slcsa/mahalanobis_file.dat c_bnnn[*]
#dump d1 all custom ${freqdump} slcsa_demo.dump id x y z c_slcsa[*]
# for testing only. in production use dump as shown above
compute max_slcsa all reduce max c_slcsa[*]
compute min_slcsa all reduce min c_slcsa[*]
thermo_style custom step pe ke temp c_max_slcsa[*] c_min_slcsa[*]
#fix 1 all nvt temp ${trequis} ${trequis} 0.100
fix 1 all npt temp ${trequis} ${trequis} 0.100 tri ${prequis_low} ${prequis_low} 1.0
run ${equilSteps}

180
examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.1 Normal file
View File

@ -0,0 +1,180 @@
LAMMPS (21 Nov 2023)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
variable trequis equal 750.0
variable prequis_low equal 0.0
variable prequis_high equal 25.0e4
variable equilSteps equal 200
variable runSteps equal 2000
variable freqdump equal 200
variable pstime equal step*dt
variable sxx equal 1.e-4*pxx
variable syy equal 1.e-4*pyy
variable szz equal 1.e-4*pzz
variable sxy equal 1.e-4*pxy
variable sxz equal 1.e-4*pxz
variable syz equal 1.e-4*pyz
variable TK equal temp
variable PE equal pe
variable KE equal ke
variable V equal vol
dimension 3
boundary p p p
units metal
atom_style atomic
read_data data.zr_cell
Reading data file ...
orthogonal box = (0 0 0) to (19.374 33.556752 30.846)
1 by 1 by 1 MPI processor grid
reading atoms ...
864 atoms
read_data CPU = 0.002 seconds
replicate 1 5 5
Replication is creating a 1x5x5 = 25 times larger system...
orthogonal box = (0 0 0) to (19.374 167.78376 154.23)
1 by 1 by 1 MPI processor grid
21600 atoms
replicate CPU = 0.001 seconds
change_box all triclinic
Changing box ...
triclinic box = (0 0 0) to (19.374 167.78376 154.23) with tilt (0 0 0)
pair_style hybrid/overlay zero 9.0 eam/fs
pair_coeff * * zero
pair_coeff * * eam/fs Zr_mm.eam.fs Zr
Reading eam/fs potential file Zr_mm.eam.fs with DATE: 2007-06-11
timestep 0.002
thermo 50
thermo_style custom step pe ke temp vol pxx pyy pzz pxy pyz pxz
# fix extra all print 50 "${pstime} ${TK} ${PE} ${KE} ${V} ${sxx} ${syy} ${szz} ${sxy} ${sxz} ${syz}" file thermo_global_npt_low_temperature_Zr_hcp.dat
velocity all create ${trequis} 42345 dist gaussian
velocity all create 750 42345 dist gaussian
# 1st step : compute the bispectrum on 24 nearest neighbors
compute bnnn all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1 delta 0.25
# 2nd step : perform dimension reduction + logistic regression
compute slcsa all slcsa/atom 8 4 dir.slcsa/mean_descriptor.dat dir.slcsa/lda_scalings.dat dir.slcsa/lr_decision.dat dir.slcsa/lr_bias.dat dir.slcsa/mahalanobis_file.dat c_bnnn[*]
Files used:
database mean descriptor: dir.slcsa/mean_descriptor.dat
lda scalings : dir.slcsa/lda_scalings.dat
lr decision : dir.slcsa/lr_decision.dat
lr bias : dir.slcsa/lr_bias.dat
maha stats : dir.slcsa/mahalanobis_file.dat
For class 0 maha threshold = 5.054
mean B:
-23.8329
4.6638
3.9805
icov:
1.1377 0.1077 -0.0171
0.1077 0.8846 -0.2577
-0.0171 -0.2577 0.6783
For class 1 maha threshold = 5.234
mean B:
-21.2853
-6.1583
1.7948
icov:
1.7124 0.0341 0.1966
0.0341 0.6453 0.288
0.1966 0.288 1.8991
For class 2 maha threshold = 5.036
mean B:
-23.1593
1.3059
-5.7549
icov:
0.7496 -0.0806 -0.1101
-0.0806 1.1178 0.1667
-0.1101 0.1667 0.6711
For class 3 maha threshold = 7.994
mean B:
68.1971
0.1604
-0.0067
icov:
0.9663 -0.1846 0.6622
-0.1846 8.2371 0.9841
0.6622 0.9841 5.9601
#dump d1 all custom ${freqdump} slcsa_demo.dump id x y z c_slcsa[*]
# for testing only. in production use dump as shown above
compute max_slcsa all reduce max c_slcsa[*]
compute min_slcsa all reduce min c_slcsa[*]
thermo_style custom step pe ke temp c_max_slcsa[*] c_min_slcsa[*]
#fix 1 all nvt temp ${trequis} ${trequis} 0.100
fix 1 all npt temp ${trequis} ${trequis} 0.100 tri ${prequis_low} ${prequis_low} 1.0
fix 1 all npt temp 750 ${trequis} 0.100 tri ${prequis_low} ${prequis_low} 1.0
fix 1 all npt temp 750 750 0.100 tri ${prequis_low} ${prequis_low} 1.0
fix 1 all npt temp 750 750 0.100 tri 0 ${prequis_low} 1.0
fix 1 all npt temp 750 750 0.100 tri 0 0 1.0
run ${equilSteps}
run 200
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11
ghost atom cutoff = 11
binsize = 5.5, bins = 4 31 29
3 neighbor lists, perpetual/occasional/extra = 2 1 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
(2) pair eam/fs, perpetual, trim from (1)
attributes: half, newton on, cut 9.6
pair build: trim
stencil: none
bin: none
(3) compute sna/atom, occasional
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 31.9 | 31.9 | 31.9 Mbytes
Step PotEng KinEng Temp c_max_slcsa[1] c_max_slcsa[2] c_max_slcsa[3] c_max_slcsa[4] c_max_slcsa[5] c_min_slcsa[1] c_min_slcsa[2] c_min_slcsa[3] c_min_slcsa[4] c_min_slcsa[5]
0 -143297.23 2093.9174 750 7.6195146 15.787294 1.2169942 111.01919 2 7.6195146 15.787294 1.2169942 111.01919 2
50 -142154.08 1007.7164 360.9442 8.8091564 19.23244 4.2093382 113.87959 2 5.0327148 9.6817454 0.02610585 106.71863 2
100 -142365.33 1406.6559 503.83647 8.6272189 17.908949 2.9294666 113.75167 2 6.2058895 11.913521 0.033775944 108.66893 2
150 -142188.18 1432.0075 512.91691 8.6441961 18.176321 2.9277374 114.27958 2 5.5899425 10.521867 0.014919473 108.14526 2
200 -142000.4 1481.7247 530.72462 8.5895692 18.65646 3.1725758 114.55015 2 5.5955774 10.776385 0.061469343 108.35384 2
Loop time of 36.3759 on 1 procs for 200 steps with 21600 atoms
Performance: 0.950 ns/day, 25.261 hours/ns, 5.498 timesteps/s, 118.760 katom-step/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 9.0837 | 9.0837 | 9.0837 | 0.0 | 24.97
Neigh | 0.52896 | 0.52896 | 0.52896 | 0.0 | 1.45
Comm | 0.045416 | 0.045416 | 0.045416 | 0.0 | 0.12
Output | 26.548 | 26.548 | 26.548 | 0.0 | 72.98
Modify | 0.1493 | 0.1493 | 0.1493 | 0.0 | 0.41
Other | | 0.02088 | | | 0.06
Nlocal: 21600 ave 21600 max 21600 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 36674 ave 36674 max 36674 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 2.61729e+06 ave 2.61729e+06 max 2.61729e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 5.24007e+06 ave 5.24007e+06 max 5.24007e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 5240069
Ave neighs/atom = 242.59579
Neighbor list builds = 4
Dangerous builds = 0
Total wall time: 0:00:43

180
examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.4 Normal file
View File

@ -0,0 +1,180 @@
LAMMPS (21 Nov 2023)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
variable trequis equal 750.0
variable prequis_low equal 0.0
variable prequis_high equal 25.0e4
variable equilSteps equal 200
variable runSteps equal 2000
variable freqdump equal 200
variable pstime equal step*dt
variable sxx equal 1.e-4*pxx
variable syy equal 1.e-4*pyy
variable szz equal 1.e-4*pzz
variable sxy equal 1.e-4*pxy
variable sxz equal 1.e-4*pxz
variable syz equal 1.e-4*pyz
variable TK equal temp
variable PE equal pe
variable KE equal ke
variable V equal vol
dimension 3
boundary p p p
units metal
atom_style atomic
read_data data.zr_cell
Reading data file ...
orthogonal box = (0 0 0) to (19.374 33.556752 30.846)
1 by 2 by 2 MPI processor grid
reading atoms ...
864 atoms
read_data CPU = 0.002 seconds
replicate 1 5 5
Replication is creating a 1x5x5 = 25 times larger system...
orthogonal box = (0 0 0) to (19.374 167.78376 154.23)
1 by 2 by 2 MPI processor grid
21600 atoms
replicate CPU = 0.001 seconds
change_box all triclinic
Changing box ...
triclinic box = (0 0 0) to (19.374 167.78376 154.23) with tilt (0 0 0)
pair_style hybrid/overlay zero 9.0 eam/fs
pair_coeff * * zero
pair_coeff * * eam/fs Zr_mm.eam.fs Zr
Reading eam/fs potential file Zr_mm.eam.fs with DATE: 2007-06-11
timestep 0.002
thermo 50
thermo_style custom step pe ke temp vol pxx pyy pzz pxy pyz pxz
# fix extra all print 50 "${pstime} ${TK} ${PE} ${KE} ${V} ${sxx} ${syy} ${szz} ${sxy} ${sxz} ${syz}" file thermo_global_npt_low_temperature_Zr_hcp.dat
velocity all create ${trequis} 42345 dist gaussian
velocity all create 750 42345 dist gaussian
# 1st step : compute the bispectrum on 24 nearest neighbors
compute bnnn all sna/atom 9.0 0.99363 8 0.5 1.0 rmin0 0.0 nnn 24 wmode 1 delta 0.25
# 2nd step : perform dimension reduction + logistic regression
compute slcsa all slcsa/atom 8 4 dir.slcsa/mean_descriptor.dat dir.slcsa/lda_scalings.dat dir.slcsa/lr_decision.dat dir.slcsa/lr_bias.dat dir.slcsa/mahalanobis_file.dat c_bnnn[*]
Files used:
database mean descriptor: dir.slcsa/mean_descriptor.dat
lda scalings : dir.slcsa/lda_scalings.dat
lr decision : dir.slcsa/lr_decision.dat
lr bias : dir.slcsa/lr_bias.dat
maha stats : dir.slcsa/mahalanobis_file.dat
For class 0 maha threshold = 5.054
mean B:
-23.8329
4.6638
3.9805
icov:
1.1377 0.1077 -0.0171
0.1077 0.8846 -0.2577
-0.0171 -0.2577 0.6783
For class 1 maha threshold = 5.234
mean B:
-21.2853
-6.1583
1.7948
icov:
1.7124 0.0341 0.1966
0.0341 0.6453 0.288
0.1966 0.288 1.8991
For class 2 maha threshold = 5.036
mean B:
-23.1593
1.3059
-5.7549
icov:
0.7496 -0.0806 -0.1101
-0.0806 1.1178 0.1667
-0.1101 0.1667 0.6711
For class 3 maha threshold = 7.994
mean B:
68.1971
0.1604
-0.0067
icov:
0.9663 -0.1846 0.6622
-0.1846 8.2371 0.9841
0.6622 0.9841 5.9601
#dump d1 all custom ${freqdump} slcsa_demo.dump id x y z c_slcsa[*]
# for testing only. in production use dump as shown above
compute max_slcsa all reduce max c_slcsa[*]
compute min_slcsa all reduce min c_slcsa[*]
thermo_style custom step pe ke temp c_max_slcsa[*] c_min_slcsa[*]
#fix 1 all nvt temp ${trequis} ${trequis} 0.100
fix 1 all npt temp ${trequis} ${trequis} 0.100 tri ${prequis_low} ${prequis_low} 1.0
fix 1 all npt temp 750 ${trequis} 0.100 tri ${prequis_low} ${prequis_low} 1.0
fix 1 all npt temp 750 750 0.100 tri ${prequis_low} ${prequis_low} 1.0
fix 1 all npt temp 750 750 0.100 tri 0 ${prequis_low} 1.0
fix 1 all npt temp 750 750 0.100 tri 0 0 1.0
run ${equilSteps}
run 200
Neighbor list info ...
update: every = 1 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 11
ghost atom cutoff = 11
binsize = 5.5, bins = 4 31 29
3 neighbor lists, perpetual/occasional/extra = 2 1 0
(1) pair zero, perpetual
attributes: half, newton on
pair build: half/bin/newton/tri
stencil: half/bin/3d/tri
bin: standard
(2) pair eam/fs, perpetual, trim from (1)
attributes: half, newton on, cut 9.6
pair build: trim
stencil: none
bin: none
(3) compute sna/atom, occasional
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 11.26 | 11.45 | 11.64 Mbytes
Step PotEng KinEng Temp c_max_slcsa[1] c_max_slcsa[2] c_max_slcsa[3] c_max_slcsa[4] c_max_slcsa[5] c_min_slcsa[1] c_min_slcsa[2] c_min_slcsa[3] c_min_slcsa[4] c_min_slcsa[5]
0 -143297.23 2093.9174 750 7.6195146 15.787294 1.2169942 111.01919 2 7.6195146 15.787294 1.2169942 111.01919 2
50 -142154.08 1007.7164 360.9442 8.8091564 19.23244 4.2093382 113.87959 2 5.0327148 9.6817454 0.02610585 106.71863 2
100 -142365.33 1406.6559 503.83647 8.6272189 17.908949 2.9294666 113.75167 2 6.2058895 11.913521 0.033775944 108.66893 2
150 -142188.18 1432.0075 512.91691 8.6441961 18.176321 2.9277374 114.27958 2 5.5899425 10.521867 0.014919473 108.14526 2
200 -142000.4 1481.7247 530.72462 8.5895692 18.65646 3.1725758 114.55015 2 5.5955774 10.776385 0.061469343 108.35384 2
Loop time of 9.81677 on 4 procs for 200 steps with 21600 atoms
Performance: 3.521 ns/day, 6.817 hours/ns, 20.373 timesteps/s, 440.063 katom-step/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.6508 | 2.6589 | 2.6698 | 0.5 | 27.09
Neigh | 0.1516 | 0.15276 | 0.15406 | 0.3 | 1.56
Comm | 0.047132 | 0.058969 | 0.066095 | 3.2 | 0.60
Output | 6.8886 | 6.8886 | 6.8886 | 0.0 | 70.17
Modify | 0.046437 | 0.04661 | 0.046825 | 0.1 | 0.47
Other | | 0.01091 | | | 0.11
Nlocal: 5400 ave 5416 max 5393 min
Histogram: 2 1 0 0 0 0 0 0 0 1
Nghost: 12902.8 ave 12911 max 12888 min
Histogram: 1 0 0 0 0 0 1 0 0 2
Neighs: 654322 ave 655602 max 650912 min
Histogram: 1 0 0 0 0 0 0 0 0 3
FullNghs: 1.31002e+06 ave 1.31507e+06 max 1.30683e+06 min
Histogram: 1 1 0 1 0 0 0 0 0 1
Total # of neighbors = 5240065
Ave neighs/atom = 242.5956
Neighbor list builds = 4
Dangerous builds = 0
Total wall time: 0:00:11

2
src/.gitignore vendored
View File

@ -635,6 +635,8 @@
/compute_rattlers_atom.h
/compute_rigid_local.cpp
/compute_rigid_local.h
/compute_slcsa_atom.cpp
/compute_slcsa_atom.h
/compute_smd_triangle_vertices.cpp
/compute_smd_triangle_vertices.h
/compute_spec_atom.cpp

416
src/EXTRA-COMPUTE/compute_slcsa_atom.cpp Normal file
View File

@ -0,0 +1,416 @@
/* ----------------------------------------------------------------------
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: Paul Lafourcade (CEA-DAM-DIF, Arpajon, France)
------------------------------------------------------------------------- */
#include "compute_slcsa_atom.h"
#include "arg_info.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "pair.h"
#include "potential_file_reader.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <iostream>
using namespace LAMMPS_NS;
static const char cite_compute_slcsa_atom_c[] =
"compute slcsa/atom command: doi:10.1088/0965-0393/21/5/055020\n\n"
"@Article{Lafourcade2023,\n"
" author = {P. Lafourcade and J.-B. Maillet and C. Denoual and E. Duval and A. Allera and A. "
"M. Goryaeva and M.-C. Marinica},\n"
" title = {Robust crystal structure identification at extreme conditions using a "
"density-independent spectral descriptor and supervised learning},\n"
" journal = {Computational Materials Science},\n"
" year = 2023,\n"
" volume = XX,\n"
" pages = {XXXXXX}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
ComputeSLCSAAtom::ComputeSLCSAAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), list(nullptr), lda_scalings(nullptr),
database_mean_descriptor(nullptr), lr_bias(nullptr), lr_decision(nullptr), icov_list(nullptr),
mean_projected_descriptors(nullptr), maha_thresholds(nullptr), full_descriptor(nullptr),
projected_descriptor(nullptr), scores(nullptr), probas(nullptr), prodright(nullptr),
dmaha(nullptr), classification(nullptr)
{
// command : compute c1 all slcsa/atom jmax nclasses parameters_file.dat
// example : compute c1 all slcsa/atom 8 4 slcsa_parameters.dat
// example : compute c1 all slcsa/atom 8 4 database_mean_descriptor.dat lda_scalings.dat lr_decision.dat lr_bias.dat mahalanobis_data.dat c_b1[*]
// Steps :
// 1. bs=bs-xbar
// 2. dred=coefs_lda*bs
// 3. scores=decision_lr*dred + lr_bias
// 4. probas=exp(scores)/sum(exp(scores))
// 5. cs=argmax(probas)
// Read the parameters file in one bloc
// File structure :
// # database mean descriptor
// vector with bso4dim rows x 1 col
// # LDA dimension reduction matrix
// matrix bso4dim rows x nclasses-1 cols
// # LR decision matrix
// matrix with nclasses rows x nclasses-1 cols
// # LR bias vector
// vector with 1 row x nclasses cols
if (narg != 11) utils::missing_cmd_args(FLERR, "compute slcsa/atom", error);
int twojmax = utils::inumeric(FLERR, arg[3], false, lmp);
if (twojmax < 0)
error->all(FLERR, "Illegal compute slcsa/atom command: twojmax must be a non-negative integer");
ncomps = compute_ncomps(twojmax);
nclasses = utils::inumeric(FLERR, arg[4], false, lmp);
if (nclasses < 2)
error->all(FLERR, "Illegal compute slcsa/atom command: nclasses must be greater than 1");
database_mean_descriptor_file = arg[5];
lda_scalings_file = arg[6];
lr_decision_file = arg[7];
lr_bias_file = arg[8];
maha_file = arg[9];
if (comm->me == 0) {
auto mesg = fmt::format(
"Files used:\n {:24}: {}\n {:24}: {}\n {:24}: {}\n {:24}: {}\n {:24}: {}\n",
"database mean descriptor", database_mean_descriptor_file, "lda scalings",
lda_scalings_file, "lr decision", lr_decision_file, "lr bias", lr_bias_file, "maha stats",
maha_file);
utils::logmesg(lmp, mesg);
}
int expand = 0;
char **earg;
int nvalues = utils::expand_args(FLERR, narg - 10, &arg[10], 1, earg, lmp);
if (earg != &arg[10]) expand = 1;
arg = earg;
ArgInfo argi(arg[0]);
value_t val;
val.id = "";
val.val.c = nullptr;
val.which = argi.get_type();
val.argindex = argi.get_index1();
val.id = argi.get_name();
if ((val.which == ArgInfo::FIX) || (val.which == ArgInfo::VARIABLE) ||
(val.which == ArgInfo::UNKNOWN) || (val.which == ArgInfo::NONE) || (argi.get_dim() > 1))
error->all(FLERR, "Invalid compute slcsa/atom argument: {}", arg[0]);
// if wildcard expansion occurred, free earg memory from exapnd_args()
if (expand) {
for (int i = 0; i < nvalues; i++) delete[] earg[i];
memory->sfree(earg);
}
val.val.c = modify->get_compute_by_id(val.id);
if (!val.val.c) error->all(FLERR, "Compute ID {} for fix slcsa/atom does not exist", val.id);
if (val.val.c->peratom_flag == 0)
error->all(FLERR, "Compute slcsa/atom compute {} does not calculate per-atom values", val.id);
if (val.argindex == 0 && val.val.c->size_peratom_cols != 0)
error->all(FLERR, "Compute slcsa/atom compute {} does not calculate a per-atom vector", val.id);
if (val.argindex && val.val.c->size_peratom_cols == 0)
error->all(FLERR, "Compute slcsa/atom compute {} does not calculate a per-atom array", val.id);
if (val.argindex && val.argindex > val.val.c->size_peratom_cols)
error->all(FLERR, "Compute slcsa/atom compute {} array is accessed out-of-range", val.id);
descriptorval = val;
memory->create(database_mean_descriptor, ncomps, "slcsa/atom:database_mean_descriptor");
memory->create(lda_scalings, ncomps, nclasses - 1, "slcsa/atom:lda_scalings");
memory->create(lr_decision, nclasses, nclasses - 1, "slcsa/atom:lr_decision");
memory->create(lr_bias, nclasses, "slcsa/atom:lr_bias");
memory->create(maha_thresholds, nclasses, "slcsa/atom:maha_thresholds");
memory->create(icov_list, nclasses, nclasses - 1, nclasses - 1, "slcsa/atom:icov_list");
memory->create(mean_projected_descriptors, nclasses, nclasses - 1,
"slcsa/atom:mean_projected_descriptors");
if (comm->me == 0) {
if (strcmp(database_mean_descriptor_file, "NULL") == 0) {
error->one(FLERR,
"Cannot open database mean descriptor file {}: ", database_mean_descriptor_file,
utils::getsyserror());
} else {
PotentialFileReader reader(lmp, database_mean_descriptor_file,
"database mean descriptor file");
int nread = 0;
while (nread < ncomps) {
auto values = reader.next_values(0);
database_mean_descriptor[nread] = values.next_double();
nread++;
}
}
if (strcmp(lda_scalings_file, "NULL") == 0) {
error->one(FLERR, "Cannot open database linear discriminant analysis scalings file {}: ",
lda_scalings_file, utils::getsyserror());
} else {
PotentialFileReader reader(lmp, lda_scalings_file, "lda scalings file");
int nread = 0;
while (nread < ncomps) {
auto values = reader.next_values(nclasses - 1);
lda_scalings[nread][0] = values.next_double();
lda_scalings[nread][1] = values.next_double();
lda_scalings[nread][2] = values.next_double();
nread++;
}
}
if (strcmp(lr_decision_file, "NULL") == 0) {
error->one(FLERR, "Cannot open logistic regression decision file {}: ", lr_decision_file,
utils::getsyserror());
} else {
PotentialFileReader reader(lmp, lr_decision_file, "lr decision file");
int nread = 0;
while (nread < nclasses) {
auto values = reader.next_values(nclasses - 1);
lr_decision[nread][0] = values.next_double();
lr_decision[nread][1] = values.next_double();
lr_decision[nread][2] = values.next_double();
nread++;
}
}
if (strcmp(lr_bias_file, "NULL") == 0) {
error->one(FLERR, "Cannot open logistic regression bias file {}: ", lr_bias_file,
utils::getsyserror());
} else {
PotentialFileReader reader(lmp, lr_bias_file, "lr bias file");
auto values = reader.next_values(nclasses);
lr_bias[0] = values.next_double();
lr_bias[1] = values.next_double();
lr_bias[2] = values.next_double();
lr_bias[3] = values.next_double();
}
if (strcmp(maha_file, "NULL") == 0) {
error->one(FLERR, "Cannot open mahalanobis stats file {}: ", maha_file, utils::getsyserror());
} else {
PotentialFileReader reader(lmp, maha_file, "mahalanobis stats file");
int nvalues = nclasses * ((nclasses - 1) * (nclasses - 1) + nclasses);
auto values = reader.next_values(nvalues);
for (int i = 0; i < nclasses; i++) {
maha_thresholds[i] = values.next_double();
for (int j = 0; j < nclasses - 1; j++)
mean_projected_descriptors[i][j] = values.next_double();
for (int k = 0; k < nclasses - 1; k++)
for (int l = 0; l < nclasses - 1; l++) icov_list[i][k][l] = values.next_double();
}
for (int i = 0; i < nclasses; i++) {
auto mesg = fmt::format("For class {} maha threshold = {:.6}\n", i, maha_thresholds[i]);
mesg += " mean B:\n";
for (int j = 0; j < nclasses - 1; j++)
mesg += fmt::format(" {:11.6}\n", mean_projected_descriptors[i][j]);
mesg += " icov:\n";
for (int j = 0; j < nclasses - 1; j++) {
mesg += fmt::format(" {:11.6} {:11.6} {:11.6}\n", icov_list[i][j][0],
icov_list[i][j][1], icov_list[i][j][2]);
}
utils::logmesg(lmp, mesg);
}
}
}
MPI_Bcast(&database_mean_descriptor[0], ncomps, MPI_DOUBLE, 0, world);
MPI_Bcast(&lda_scalings[0][0], ncomps * (nclasses - 1), MPI_DOUBLE, 0, world);
MPI_Bcast(&lr_decision[0][0], nclasses * (nclasses - 1), MPI_DOUBLE, 0, world);
MPI_Bcast(&lr_bias[0], nclasses, MPI_DOUBLE, 0, world);
MPI_Bcast(&maha_thresholds[0], nclasses, MPI_DOUBLE, 0, world);
MPI_Bcast(&mean_projected_descriptors[0][0], nclasses * (nclasses - 1), MPI_DOUBLE, 0, world);
MPI_Bcast(&icov_list[0][0][0], nclasses * (nclasses - 1) * (nclasses - 1), MPI_DOUBLE, 0, world);
peratom_flag = 1;
size_peratom_cols = nclasses + 1;
ncols = nclasses + 1;
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeSLCSAAtom::~ComputeSLCSAAtom()
{
memory->destroy(classification);
memory->destroy(database_mean_descriptor);
memory->destroy(lda_scalings);
memory->destroy(lr_decision);
memory->destroy(lr_bias);
memory->destroy(maha_thresholds);
memory->destroy(mean_projected_descriptors);
memory->destroy(icov_list);
memory->destroy(full_descriptor);
memory->destroy(projected_descriptor);
memory->destroy(scores);
memory->destroy(probas);
memory->destroy(prodright);
memory->destroy(dmaha);
}
/* ---------------------------------------------------------------------- */
void ComputeSLCSAAtom::init()
{
if (modify->get_compute_by_style(style).size() > 1)
if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style);
}
/* ---------------------------------------------------------------------- */
void ComputeSLCSAAtom::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeSLCSAAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
// grow per-atom if necessary
if (atom->nmax > nmax) {
memory->destroy(classification);
nmax = atom->nmax;
memory->create(classification, nmax, ncols, "slcsa/atom:classification");
array_atom = classification;
}
int *mask = atom->mask;
int nlocal = atom->nlocal;
double **compute_array;
if (descriptorval.which == ArgInfo::COMPUTE) {
if (!(descriptorval.val.c->invoked_flag & Compute::INVOKED_PERATOM)) {
descriptorval.val.c->compute_peratom();
descriptorval.val.c->invoked_flag |= Compute::INVOKED_PERATOM;
}
compute_array = descriptorval.val.c->array_atom;
}
memory->create(full_descriptor, ncomps, "slcsa/atom:local descriptor");
memory->create(projected_descriptor, nclasses - 1, "slcsa/atom:reduced descriptor");
memory->create(scores, nclasses, "slcsa/atom:scores");
memory->create(probas, nclasses, "slcsa/atom:probas");
memory->create(prodright, nclasses - 1, "slcsa/atom:prodright");
memory->create(dmaha, nclasses, "slcsa/atom:prodright");
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
for (int j = 0; j < ncomps; j++) { full_descriptor[j] = compute_array[i][j]; }
// Here comes the LDA + LR process
// 1st step : Retrieve mean database descriptor
for (int j = 0; j < ncomps; j++) { full_descriptor[j] -= database_mean_descriptor[j]; }
// 2nd step : Matrix multiplication to go from ncompsx1 -> (nclasses-1)*1
for (int j = 0; j < nclasses - 1; j++) {
projected_descriptor[j] = 0.;
for (int k = 0; k < ncomps; k++) {
projected_descriptor[j] += full_descriptor[k] * lda_scalings[k][j];
}
}
// 3rd step : Matrix multiplication
for (int j = 0; j < nclasses; j++) {
scores[j] = lr_bias[j];
for (int k = 0; k < nclasses - 1; k++) {
scores[j] += lr_decision[j][k] * projected_descriptor[k];
}
}
// 4th step : Matrix multiplication
double sumexpscores = 0.;
for (int j = 0; j < nclasses; j++) sumexpscores += exp(scores[j]);
for (int j = 0; j < nclasses; j++) { probas[j] = exp(scores[j]) / sumexpscores; }
classification[i][nclasses] = argmax(probas, nclasses);
// 5th step : Mahalanobis distance
for (int j = 0; j < nclasses; j++) {
prodright[0] = 0.;
prodright[1] = 0.;
prodright[2] = 0.;
for (int k = 0; k < nclasses - 1; k++) {
for (int l = 0; l < nclasses - 1; l++) {
prodright[k] +=
(icov_list[j][k][l] * (projected_descriptor[k] - mean_projected_descriptors[j][k]));
}
}
double prodleft = 0.;
for (int k = 0; k < nclasses - 1; k++) {
prodleft += (prodright[k] * (projected_descriptor[k] - mean_projected_descriptors[j][k]));
}
classification[i][j] = sqrt(prodleft);
}
// 6th step : Sanity check
int locclass = classification[i][nclasses];
if (classification[i][locclass] > maha_thresholds[locclass]) {
classification[i][nclasses] = -1.;
}
} else {
for (int j = 0; j < ncols; j++) { classification[i][j] = -1.; }
}
}
memory->destroy(full_descriptor);
memory->destroy(projected_descriptor);
memory->destroy(scores);
memory->destroy(probas);
memory->destroy(prodright);
memory->destroy(dmaha);
}
int ComputeSLCSAAtom::compute_ncomps(int twojmax)
{
int ncount;
ncount = 0;
for (int j1 = 0; j1 <= twojmax; j1++)
for (int j2 = 0; j2 <= j1; j2++)
for (int j = j1 - j2; j <= MIN(twojmax, j1 + j2); j += 2)
if (j >= j1) ncount++;
return ncount;
}
int ComputeSLCSAAtom::argmax(double arr[], int size)
{
int maxIndex = 0; // Initialize the index of the maximum value to the first element.
double maxValue = arr[0]; // Initialize the maximum value to the first element.
for (int i = 1; i < size; ++i) {
if (arr[i] > maxValue) {
// If a greater value is found, update the maxIndex and maxValue.
maxIndex = i;
maxValue = arr[i];
}
}
return maxIndex;
}

95
src/EXTRA-COMPUTE/compute_slcsa_atom.h Normal file
View File

@ -0,0 +1,95 @@
/* -*- 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: Paul Lafourcade (CEA-DAM-DIF, Arpajon, France)
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(slcsa/atom,ComputeSLCSAAtom);
// clang-format on
#else
#ifndef LMP_COMPUTE_SLCSA_ATOM_H
#define LMP_COMPUTE_SLCSA_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeSLCSAAtom : public Compute {
public:
ComputeSLCSAAtom(class LAMMPS *, int, char **);
~ComputeSLCSAAtom() override;
void init() override;
void init_list(int, class NeighList *) override;
void compute_peratom() override;
// double memory_usage() override;
int compute_ncomps(int);
int argmax(double *, int);
private:
struct value_t {
int which; // type of data: COMPUTE, FIX, VARIABLE
int argindex; // 1-based index if data is vector, else 0
std::string id; // compute/fix/variable ID
union {
class Compute *c;
class Fix *f;
int v;
} val;
};
value_t descriptorval;
int nmax;
int ncols;
int nevery;
int ncomps;
int nclasses;
const char *database_mean_descriptor_file;
const char *lda_scalings_file;
const char *lr_decision_file;
const char *lr_bias_file;
const char *covmat_file;
const char *maha_file;
class NeighList *list;
// LDA dimension reduction
double **lda_scalings;
double *database_mean_descriptor;
// LR classification
double *lr_bias;
double **lr_decision;
// Mahalanobis distance calculation
double ***icov_list;
double **mean_projected_descriptors;
double *maha_thresholds;
// Per-atom local arrays
double *full_descriptor;
double *projected_descriptor;
double *scores;
double *probas;
double *prodright;
double *dmaha;
// Output array
double **classification;
};
} // namespace LAMMPS_NS
#endif
#endif

427
src/ML-SNAP/compute_sna_atom.cpp
View File

@ -56,7 +56,10 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
wselfallflag = 0;
switchinnerflag = 0;
nelements = 1;
nnn = 12;
wmode = 0;
delta = 1.e-3;
nearest_neighbors_mode = false;
// process required arguments
memory->create(radelem, ntypes + 1, "sna/atom:radelem"); // offset by 1 to match up with types
@ -114,6 +117,22 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
quadraticflag = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
iarg += 2;
} else if (strcmp(arg[iarg],"nnn") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
nnn = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
nearest_neighbors_mode = true;
if (nnn <= 0) error->all(FLERR, "Illegal compute compute {} command", style);
iarg += 2;
} else if (strcmp(arg[iarg],"wmode") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
wmode = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
if (wmode < 0) error->all(FLERR, "Illegal compute compute {} command", style);
iarg += 2;
} else if (strcmp(arg[iarg],"delta") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
delta = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
if (delta < 1.0e-3) error->all(FLERR, "Illegal compute compute {} command", style);
iarg += 2;
} else if (strcmp(arg[iarg], "chem") == 0) {
if (iarg + 2 > narg) error->all(FLERR, "Illegal compute {} command", style);
chemflag = 1;
@ -183,6 +202,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
nmax = 0;
sna = nullptr;
}
/* ---------------------------------------------------------------------- */
@ -209,7 +229,7 @@ void ComputeSNAAtom::init()
{
if (force->pair == nullptr)
error->all(FLERR,"Compute sna/atom requires a pair style be defined");
rcutsq = force->pair->cutforce * force->pair->cutforce;
if (cutmax > force->pair->cutforce)
error->all(FLERR,"Compute sna/atom cutoff is longer than pairwise cutoff");
@ -275,63 +295,164 @@ void ComputeSNAAtom::compute_peratom()
const int* const jlist = firstneigh[i];
const int jnum = numneigh[i];
// ensure rij, inside, and typej are of size jnum
snaptr->grow_rij(jnum);
// ############################################################################## //
// ##### Start of section for computing bispectrum on nnn nearest neighbors ##### //
// ############################################################################## //
if (nearest_neighbors_mode == true) {
// ##### 1) : consider full neighbor list in rlist
memory->create(distsq, jnum, "snann/atom:distsq");
memory->create(rlist, jnum, 3, "snann/atom:rlist");
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// typej = types of neighbors of I within cutoff
int ncount = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
int jtype = type[j];
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx * delx + dely * dely + delz * delz;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
int jelem = 0;
if (chemflag)
jelem = map[jtype];
if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jtype];
snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
if (switchinnerflag) {
snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]);
snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]);
if (rsq < rcutsq) {
distsq[ncount] = rsq;
rlist[ncount][0] = delx;
rlist[ncount][1] = dely;
rlist[ncount][2] = delz;
ncount++;
}
if (chemflag) snaptr->element[ninside] = jelem;
ninside++;
}
// ##### 2) : compute optimal cutoff such that sum weights S_target = nnn
double S_target=1.*nnn;
double rc_start=0.1;
double rc_max=sqrt(rcutsq);
double tol=1.e-8;
double * sol_dich = dichotomie(S_target, rc_start, rc_max, tol, distsq, ncount, wmode, delta);
memory->destroy(distsq);
// ##### 3) : assign that optimal cutoff radius to bispectrum context using rcsol
double rcsol = (sol_dich[0]+sol_dich[1])/2.;
memory->destroy(sol_dich);
snaptr->grow_rij(ncount);
int ninside = 0;
for (int jj = 0; jj < ncount; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double rsq = rlist[jj][0]*rlist[jj][0]+rlist[jj][1]*rlist[jj][1]+rlist[jj][2]*rlist[jj][2];
int jtype = type[j];
int jelem = 0;
if (chemflag)
jelem = map[jtype];
if (rsq < rcsol*rcsol) {
snaptr->rij[ninside][0] = rlist[jj][0];//rijmax;
snaptr->rij[ninside][1] = rlist[jj][1];//rijmax;
snaptr->rij[ninside][2] = rlist[jj][2];//rijmax;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = 1.;
snaptr->rcutij[ninside] = rcsol;
if (switchinnerflag) {
snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]);
snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]);
}
if (chemflag) snaptr->element[ninside] = jelem;
ninside++;
}
}
memory->destroy(rlist);
// ############################################################################ //
// ##### End of section for computing bispectrum on nnn nearest neighbors ##### //
// ############################################################################ //
snaptr->compute_ui(ninside, ielem);
snaptr->compute_zi();
snaptr->compute_bi(ielem);
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
sna[i][icoeff] = snaptr->blist[icoeff];
if (quadraticflag) {
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bi = snaptr->blist[icoeff];
// diagonal element of quadratic matrix
sna[i][ncount++] = 0.5*bi*bi;
// upper-triangular elements of quadratic matrix
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
sna[i][ncount++] = bi*snaptr->blist[jcoeff];
}
}
} else {
// ensure rij, inside, and typej are of size jnum
snaptr->grow_rij(jnum);
// rij[][3] = displacements between atom I and those neighbors
// inside = indices of neighbors of I within cutoff
// typej = types of neighbors of I within cutoff
int ninside = 0;
for (int jj = 0; jj < jnum; jj++) {
int j = jlist[jj];
j &= NEIGHMASK;
const double delx = xtmp - x[j][0];
const double dely = ytmp - x[j][1];
const double delz = ztmp - x[j][2];
const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j];
int jelem = 0;
if (chemflag)
jelem = map[jtype];
if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
snaptr->rij[ninside][0] = delx;
snaptr->rij[ninside][1] = dely;
snaptr->rij[ninside][2] = delz;
snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jtype];
snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac;
if (switchinnerflag) {
snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]);
snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]);
}
if (chemflag) snaptr->element[ninside] = jelem;
ninside++;
}
}
snaptr->compute_ui(ninside, ielem);
snaptr->compute_zi();
snaptr->compute_bi(ielem);
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
sna[i][icoeff] = snaptr->blist[icoeff];
if (quadraticflag) {
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bi = snaptr->blist[icoeff];
// diagonal element of quadratic matrix
sna[i][ncount++] = 0.5*bi*bi;
// upper-triangular elements of quadratic matrix
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
sna[i][ncount++] = bi*snaptr->blist[jcoeff];
}
}
}
snaptr->compute_ui(ninside, ielem);
snaptr->compute_zi();
snaptr->compute_bi(ielem);
for (int icoeff = 0; icoeff < ncoeff; icoeff++)
sna[i][icoeff] = snaptr->blist[icoeff];
if (quadraticflag) {
int ncount = ncoeff;
for (int icoeff = 0; icoeff < ncoeff; icoeff++) {
double bi = snaptr->blist[icoeff];
// diagonal element of quadratic matrix
sna[i][ncount++] = 0.5*bi*bi;
// upper-triangular elements of quadratic matrix
for (int jcoeff = icoeff+1; jcoeff < ncoeff; jcoeff++)
sna[i][ncount++] = bi*snaptr->blist[jcoeff];
}
}
} else {
for (int icoeff = 0; icoeff < size_peratom_cols; icoeff++)
sna[i][icoeff] = 0.0;
@ -352,3 +473,207 @@ double ComputeSNAAtom::memory_usage()
return bytes;
}
/* ----------------------------------------------------------------------
select3 routine from Numerical Recipes (slightly modified)
find k smallest values in array of length n
sort auxiliary arrays at same time
------------------------------------------------------------------------- */
// Use no-op do while to create single statement
#define SWAP(a, b) \
do { \
tmp = a; \
(a) = b; \
(b) = tmp; \
} while (0)
#define ISWAP(a, b) \
do { \
itmp = a; \
(a) = b; \
(b) = itmp; \
} while (0)
#define SWAP3(a, b) \
do { \
tmp = (a)[0]; \
(a)[0] = (b)[0]; \
(b)[0] = tmp; \
tmp = (a)[1]; \
(a)[1] = (b)[1]; \
(b)[1] = tmp; \
tmp = (a)[2]; \
(a)[2] = (b)[2]; \
(b)[2] = tmp; \
} while (0)
/* ---------------------------------------------------------------------- */
void ComputeSNAAtom::select3(int k, int n, double *arr, int *iarr, double **arr3)
{
int i, ir, j, l, mid, ia, itmp;
double a, tmp, a3[3];
arr--;
iarr--;
arr3--;
l = 1;
ir = n;
for (;;) {
if (ir <= l + 1) {
if (ir == l + 1 && arr[ir] < arr[l]) {
SWAP(arr[l], arr[ir]);
ISWAP(iarr[l], iarr[ir]);
SWAP3(arr3[l], arr3[ir]);
}
return;
} else {
mid = (l + ir) >> 1;
SWAP(arr[mid], arr[l + 1]);
ISWAP(iarr[mid], iarr[l + 1]);
SWAP3(arr3[mid], arr3[l + 1]);
if (arr[l] > arr[ir]) {
SWAP(arr[l], arr[ir]);
ISWAP(iarr[l], iarr[ir]);
SWAP3(arr3[l], arr3[ir]);
}
if (arr[l + 1] > arr[ir]) {
SWAP(arr[l + 1], arr[ir]);
ISWAP(iarr[l + 1], iarr[ir]);
SWAP3(arr3[l + 1], arr3[ir]);
}
if (arr[l] > arr[l + 1]) {
SWAP(arr[l], arr[l + 1]);
ISWAP(iarr[l], iarr[l + 1]);
SWAP3(arr3[l], arr3[l + 1]);
}
i = l + 1;
j = ir;
a = arr[l + 1];
ia = iarr[l + 1];
a3[0] = arr3[l + 1][0];
a3[1] = arr3[l + 1][1];
a3[2] = arr3[l + 1][2];
for (;;) {
do i++;
while (arr[i] < a);
do j--;
while (arr[j] > a);
if (j < i) break;
SWAP(arr[i], arr[j]);
ISWAP(iarr[i], iarr[j]);
SWAP3(arr3[i], arr3[j]);
}
arr[l + 1] = arr[j];
arr[j] = a;
iarr[l + 1] = iarr[j];
iarr[j] = ia;
arr3[l + 1][0] = arr3[j][0];
arr3[l + 1][1] = arr3[j][1];
arr3[l + 1][2] = arr3[j][2];
arr3[j][0] = a3[0];
arr3[j][1] = a3[1];
arr3[j][2] = a3[2];
if (j >= k) ir = j - 1;
if (j <= k) l = i;
}
}
}
double * ComputeSNAAtom::weights(double * rsq, double rcut, int ncounts)
{
double * w=nullptr;
memory->destroy(w);
memory->create(w, ncounts, "snann:gauss_weights");
double rloc=0.;
for (int i=0; i<ncounts; i++)
{
rloc = sqrt(rsq[i]);
if (rloc > rcut){
w[i]=0.;
} else {
w[i]=1.;
}
}
return w;
}
double * ComputeSNAAtom::tanh_weights(double * rsq, double rcut, double delta, int ncounts)
{
double * w=nullptr;
memory->destroy(w);
memory->create(w, ncounts, "snann:gauss_weights");
double rloc=0.;
for (int i=0; i<ncounts; i++)
{
rloc = sqrt(rsq[i]);
w[i] = 0.5*(1.-tanh((rloc-rcut)/delta));
}
return w;
}
double ComputeSNAAtom::sum_weights(double * rsq, double * w, int ncounts)
{
double S=0.;
double rloc=0.;
for (int i=0; i<ncounts; i++)
{
S += w[i];
}
return S;
}
double ComputeSNAAtom::get_target_rcut(double S_target, double * rsq, double rcut, int ncounts, int weightmode, double delta)
{
double S_sol;
if (weightmode == 0) {
double * www = weights(rsq, rcut, ncounts);
S_sol = sum_weights(rsq, www, ncounts);
memory->destroy(www);
} else if (weightmode == 1) {
double * www = tanh_weights(rsq, rcut, delta, ncounts);
S_sol = sum_weights(rsq, www, ncounts);
memory->destroy(www);
}
double err = S_sol - S_target;
return err;
}
double * ComputeSNAAtom::dichotomie(double S_target, double a, double b, double e, double * rsq, int ncounts, int weightmode, double delta)
{
double d=b-a;
double * sol = nullptr;
memory->destroy(sol);
memory->create(sol, 2, "snann:sol");
double m=0.;
int cnt=0;
do
{
m = ( a + b ) / 2.;
d = fabs( b - a );
double f_ra = get_target_rcut(S_target, rsq, a, ncounts, weightmode, delta);
double f_rm = get_target_rcut(S_target, rsq, m, ncounts, weightmode, delta);
if (f_rm == 0.)
{
sol[0]=m;
sol[1]=m;
return sol;
}
else if (f_rm*f_ra > 0.)
{
a = m;
}
else
{
b = m;
}
cnt+=1;
} while ( d > e );
sol[0]=a;
sol[1]=b;
return sol;
}

19
src/ML-SNAP/compute_sna_atom.h
View File

@ -11,6 +11,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Paul Lafourcade (CEA-DAM-DIF, Arpajon, France)
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(sna/atom,ComputeSNAAtom);
@ -32,10 +36,25 @@ class ComputeSNAAtom : public Compute {
void init_list(int, class NeighList *) override;
void compute_peratom() override;
double memory_usage() override;
double rcutsq;
void select3(int, int, double *, int *, double **);
double * weights(double *, double, int);
double * tanh_weights(double *, double, double, int);
double sum_weights(double *, double *, int);
double get_target_rcut(double, double *, double, int, int, double);
double * dichotomie(double, double, double, double, double *, int, int, double);
private:
int nmax;
int ncoeff;
int nnn;
int wmode;
double delta;
bool nearest_neighbors_mode;
double *distsq;
double **rlist;
int *nearest;
double **cutsq;
class NeighList *list;
double **sna;