diff --git a/doc/src/Bibliography.rst b/doc/src/Bibliography.rst index e9ea8b0925..4ed8e73dfe 100644 --- a/doc/src/Bibliography.rst +++ b/doc/src/Bibliography.rst @@ -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) `_ + **(Lamoureux and Roux)** G.\ Lamoureux, B. Roux, J. Chem. Phys 119, 3025 (2003) diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index b7e2156d88..ab4d8d64a2 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -123,6 +123,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`reduce/region ` * :doc:`rigid/local ` * :doc:`saed ` + * :doc:`slcsa/atom ` * :doc:`slice ` * :doc:`smd/contact/radius ` * :doc:`smd/damage ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6491696271..afe57d9048 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -287,6 +287,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`reduce/region ` - same as compute reduce, within a region * :doc:`rigid/local ` - extract rigid body attributes * :doc:`saed ` - electron diffraction intensity on a mesh of reciprocal lattice nodes +* :doc:`slcsa/atom ` - perform Supervised Learning Crystal Structure Analysis (SL-CSA) * :doc:`slice ` - extract values from global vector or array * :doc:`smd/contact/radius ` - contact radius for Smooth Mach Dynamics * :doc:`smd/damage ` - damage status of SPH particles in Smooth Mach Dynamics diff --git a/doc/src/compute_slcsa_atom.rst b/doc/src/compute_slcsa_atom.rst new file mode 100644 index 0000000000..6b2708c4d9 --- /dev/null +++ b/doc/src/compute_slcsa_atom.rst @@ -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 ` 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) ` +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 `_. + +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 `_ 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 ` 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 ` 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 +` page for more info. + +Related commands +"""""""""""""""" + +:doc:`compute sna/atom ` + +Default +""""""" + +none + +---------- + +.. _Lafourcade2023_1: + +**(Lafourcade)** Lafourcade, Maillet, Denoual, Duval, Allera, Goryaeva, and Marinica, +`Comp. Mat. Science, 230, 112534 (2023) `_ diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index 8d06868f3d..179c362dc6 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -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 ` 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) `. + ---------- Output info @@ -585,6 +615,7 @@ Related commands """""""""""""""" :doc:`pair_style snap ` +:doc:`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) `_ diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index c31c6a0a70..272e747358 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/Zr_mm.eam.fs b/examples/PACKAGES/sna_nnn_slcsa/Zr_mm.eam.fs new file mode 120000 index 0000000000..f8b779307d --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/Zr_mm.eam.fs @@ -0,0 +1 @@ +../../../potentials/Zr_mm.eam.fs \ No newline at end of file diff --git a/examples/PACKAGES/sna_nnn_slcsa/data.zr_cell b/examples/PACKAGES/sna_nnn_slcsa/data.zr_cell new file mode 100644 index 0000000000..d7c83fc716 --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/data.zr_cell @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lda_scalings.dat b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lda_scalings.dat new file mode 100644 index 0000000000..68a78f8c40 --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lda_scalings.dat @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_bias.dat b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_bias.dat new file mode 100644 index 0000000000..467f1844d9 --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_bias.dat @@ -0,0 +1 @@ +-6.32012657 5.62127377 1.19871662 -0.49986382 diff --git a/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_decision.dat b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_decision.dat new file mode 100644 index 0000000000..e938d59d8b --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/lr_decision.dat @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mahalanobis_file.dat b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mahalanobis_file.dat new file mode 100644 index 0000000000..299ef3c72f --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mahalanobis_file.dat @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mean_descriptor.dat b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mean_descriptor.dat new file mode 100644 index 0000000000..ae50ec809f --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/dir.slcsa/mean_descriptor.dat @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/in.slcsa b/examples/PACKAGES/sna_nnn_slcsa/in.slcsa new file mode 100644 index 0000000000..31a8189da3 --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/in.slcsa @@ -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} diff --git a/examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.1 b/examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.1 new file mode 100644 index 0000000000..58c2f40684 --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.1 @@ -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 diff --git a/examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.4 b/examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.4 new file mode 100644 index 0000000000..6436eabe2b --- /dev/null +++ b/examples/PACKAGES/sna_nnn_slcsa/log.12Dec23.slcsa.g++.4 @@ -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 diff --git a/src/.gitignore b/src/.gitignore index 74164574c2..3b7902303d 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -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 diff --git a/src/EXTRA-COMPUTE/compute_slcsa_atom.cpp b/src/EXTRA-COMPUTE/compute_slcsa_atom.cpp new file mode 100644 index 0000000000..f12551085c --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_slcsa_atom.cpp @@ -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 +#include +#include + +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; +} diff --git a/src/EXTRA-COMPUTE/compute_slcsa_atom.h b/src/EXTRA-COMPUTE/compute_slcsa_atom.h new file mode 100644 index 0000000000..6d7cd90c31 --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_slcsa_atom.h @@ -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 diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index 2de25b09b6..753751690d 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -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 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; idestroy(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; +} diff --git a/src/ML-SNAP/compute_sna_atom.h b/src/ML-SNAP/compute_sna_atom.h index 29a84c8dcf..2283865431 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -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;