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

Merge remote-tracking branch 'refs/remotes/origin/chem_snap' into chem_snap

Browse Source
This commit is contained in:
Sievers
2020-06-15 12:22:11 -07:00
parent 39aee089e7 5d3f67dff4
commit ac87f1763c
28 changed files with 626 additions and 167 deletions
Download Patch File Download Diff File Expand all files Collapse all files

154
doc/src/compute_sna_atom.rst
View File

@ -30,7 +30,7 @@ Syntax
* R_1, R_2,... = list of cutoff radii, one for each type (distance units) * R_1, R_2,... = list of cutoff radii, one for each type (distance units)
* w_1, w_2,... = list of neighbor weights, one for each type * w_1, w_2,... = list of neighbor weights, one for each type
* zero or more keyword/value pairs may be appended * zero or more keyword/value pairs may be appended
* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* * keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag*
.. parsed-literal:: .. parsed-literal::
@ -44,6 +44,15 @@ Syntax
*quadraticflag* value = *0* or *1* *quadraticflag* value = *0* or *1*
*0* = do not generate quadratic terms *0* = do not generate quadratic terms
*1* = generate quadratic terms *1* = generate quadratic terms
*chem* values = *nelements* *elementlist*
*nelements* = number of SNAP elements
*elementlist* = *ntypes* integers in range [0, *nelements*)
*bnormflag* value = *0* or *1*
*0* = do not normalize
*1* = normalize bispectrum components
*wselfallflag* value = *0* or *1*
*0* = self-contribution only for element of central atom
*1* = self-contribution for all elements
Examples Examples
"""""""" """"""""
@ -54,6 +63,7 @@ Examples
compute db all sna/atom 1.4 0.95 6 2.0 1.0 compute db all sna/atom 1.4 0.95 6 2.0 1.0
compute vb all sna/atom 1.4 0.95 6 2.0 1.0 compute vb all sna/atom 1.4 0.95 6 2.0 1.0
compute snap all snap 1.4 0.95 6 2.0 1.0 compute snap all snap 1.4 0.95 6 2.0 1.0
compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 chem 2 0 1
Description Description
""""""""""" """""""""""
@ -71,27 +81,26 @@ mathematical definition is given in the paper by Thompson et
al. :ref:`(Thompson) <Thompson20141>` al. :ref:`(Thompson) <Thompson20141>`
The position of a neighbor atom *i'* relative to a central atom *i* is The position of a neighbor atom *i'* relative to a central atom *i* is
a point within the 3D ball of radius *R_ii' = rcutfac\*(R_i + R_i')* a point within the 3D ball of radius :math:`R_{ii'}` = *rcutfac* :math:`(R_i + R_i')`
Bartok et al. :ref:`(Bartok) <Bartok20101>`, proposed mapping this 3D ball Bartok et al. :ref:`(Bartok) <Bartok20101>`, proposed mapping this 3D ball
onto the 3-sphere, the surface of the unit ball in a four-dimensional onto the 3-sphere, the surface of the unit ball in a four-dimensional
space. The radial distance *r* within *R_ii'* is mapped on to a third space. The radial distance *r* within *R_ii'* is mapped on to a third
polar angle *theta0* defined by, polar angle :math:`\theta_0` defined by,
.. math:: .. math::
\theta_0 = {\tt rfac0} \frac{r-r_{min0}}{R_{ii'}-r_{min0}} \pi \theta_0 = {\sf rfac0} \frac{r-r_{min0}}{R_{ii'}-r_{min0}} \pi
In this way, all possible neighbor positions are mapped on to a subset In this way, all possible neighbor positions are mapped on to a subset
of the 3-sphere. Points south of the latitude *theta0max=rfac0\*Pi* of the 3-sphere. Points south of the latitude :math:`\theta_0` = *rfac0* :math:`\pi`
are excluded. are excluded.
The natural basis for functions on the 3-sphere is formed by the 4D The natural basis for functions on the 3-sphere is formed by the
hyperspherical harmonics *U\^j_m,m'(theta, phi, theta0).* These representatives of *SU(2)*, the matrices :math:`U^j_{m,m'}(\theta, \phi, \theta_0)`.
functions are better known as *D\^j_m,m',* the elements of the Wigner These functions are better known as :math:`D^j_{m,m'}`, the elements of the Wigner
*D*\ -matrices :ref:`(Meremianin <Meremianin2006>`, *D*\ -matrices :ref:`(Meremianin <Meremianin2006>`,
:ref:`Varshalovich) <Varshalovich1987>`. :ref:`Varshalovich <Varshalovich1987>`, :ref:`Mason) <Mason2009>`
The density of neighbors on the 3-sphere can be written as a sum of The density of neighbors on the 3-sphere can be written as a sum of
Dirac-delta functions, one for each neighbor, weighted by species and Dirac-delta functions, one for each neighbor, weighted by species and
radial distance. Expanding this density function as a generalized radial distance. Expanding this density function as a generalized
@ -100,20 +109,20 @@ coefficient as
.. math:: .. math::
u^j_{m,m'} = U^j_{m,m'}(0,0,0) + \sum_{r_{ii'} < R_{ii'}}{f_c(r_{ii'}) w_{i'} U^j_{m,m'}(\theta_0,\theta,\phi)} u^j_{m,m'} = U^j_{m,m'}(0,0,0) + \sum_{r_{ii'} < R_{ii'}}{f_c(r_{ii'}) w_{\mu_{i'}} U^j_{m,m'}(\theta_0,\theta,\phi)}
The *w_i'* neighbor weights are dimensionless numbers that are chosen The :math:`w_{\mu_{i'}}` neighbor weights are dimensionless numbers that depend on
to distinguish atoms of different types, while the central atom is :math:`\mu_{i'}`, the SNAP element of atom *i'*, while the central atom is
arbitrarily assigned a unit weight. The function *fc(r)* ensures that arbitrarily assigned a unit weight. The function :math:`f_c(r)` ensures that
the contribution of each neighbor atom goes smoothly to zero at the contribution of each neighbor atom goes smoothly to zero at
*R_ii'*: :math:`R_{ii'}`:
.. math:: .. math::
f_c(r) = & \frac{1}{2}(\cos(\pi \frac{r-r_{min0}}{R_{ii'}-r_{min0}}) + 1), r \leq R_{ii'} \\ f_c(r) = & \frac{1}{2}(\cos(\pi \frac{r-r_{min0}}{R_{ii'}-r_{min0}}) + 1), r \leq R_{ii'} \\
= & 0, r > R_{ii'} = & 0, r > R_{ii'}
The expansion coefficients *u\^j_m,m'* are complex-valued and they are The expansion coefficients :math:`u^j_{m,m'}` are complex-valued and they are
not directly useful as descriptors, because they are not invariant not directly useful as descriptors, because they are not invariant
under rotation of the polar coordinate frame. However, the following under rotation of the polar coordinate frame. However, the following
scalar triple products of expansion coefficients can be shown to be scalar triple products of expansion coefficients can be shown to be
@ -128,7 +137,8 @@ real-valued and invariant under rotation :ref:`(Bartok) <Bartok20101>`.
{j_2} {m_2} {m'_2} \end{array}} {j_2} {m_2} {m'_2} \end{array}}
u^{j_1}_{m_1,m'_1} u^{j_2}_{m_2,m'_2} u^{j_1}_{m_1,m'_1} u^{j_2}_{m_2,m'_2}
The constants *H\^jmm'_j1m1m1'_j2m2m2'* are coupling coefficients, The constants :math:`H^{jmm'}_{j_1 m_1 m_{1'},j_2 m_ 2m_{2'}}`
are coupling coefficients,
analogous to Clebsch-Gordan coefficients for rotations on the analogous to Clebsch-Gordan coefficients for rotations on the
2-sphere. These invariants are the components of the bispectrum and 2-sphere. These invariants are the components of the bispectrum and
these are the quantities calculated by the compute *sna/atom*\ . They these are the quantities calculated by the compute *sna/atom*\ . They
@ -136,13 +146,12 @@ characterize the strength of density correlations at three points on
the 3-sphere. The j2=0 subset form the power spectrum, which the 3-sphere. The j2=0 subset form the power spectrum, which
characterizes the correlations of two points. The lowest-order characterizes the correlations of two points. The lowest-order
components describe the coarsest features of the density function, components describe the coarsest features of the density function,
while higher-order components reflect finer detail. Note that the while higher-order components reflect finer detail. Each bispectrum
central atom is included in the expansion, so three point-correlations component contains terms that depend on the positions of up to 4
can be either due to three neighbors, or two neighbors and the central atoms (3 neighbors and the central atom).
atom.
Compute *snad/atom* calculates the derivative of the bispectrum components Compute *snad/atom* calculates the derivative of the bispectrum components
summed separately for each atom type: summed separately for each LAMMPS atom type:
.. math:: .. math::
@ -165,7 +174,7 @@ Again, the sum is over all atoms *i'* of atom type *I*\ . For each atom
virial components, each atom type, and each bispectrum component. See virial components, each atom type, and each bispectrum component. See
section below on output for a detailed explanation. section below on output for a detailed explanation.
Compute *snap* calculates a global array contains information related Compute *snap* calculates a global array containing information related
to all three of the above per-atom computes *sna/atom*\ , *snad/atom*\ , to all three of the above per-atom computes *sna/atom*\ , *snad/atom*\ ,
and *snav/atom*\ . The first row of the array contains the summation of and *snav/atom*\ . The first row of the array contains the summation of
*sna/atom* over all atoms, but broken out by type. The last six rows *sna/atom* over all atoms, but broken out by type. The last six rows
@ -201,8 +210,8 @@ The argument *rcutfac* is a scale factor that controls the ratio of
atomic radius to radial cutoff distance. atomic radius to radial cutoff distance.
The argument *rfac0* and the optional keyword *rmin0* define the The argument *rfac0* and the optional keyword *rmin0* define the
linear mapping from radial distance to polar angle *theta0* on the linear mapping from radial distance to polar angle :math:`theta_0` on the
3-sphere. 3-sphere, given above.
The argument *twojmax* defines which The argument *twojmax* defines which
bispectrum components are generated. See section below on output for a bispectrum components are generated. See section below on output for a
@ -210,7 +219,7 @@ detailed explanation of the number of bispectrum components and the
ordered in which they are listed. ordered in which they are listed.
The keyword *switchflag* can be used to turn off the switching The keyword *switchflag* can be used to turn off the switching
function. function :math:`f_c(r)`.
The keyword *bzeroflag* determines whether or not *B0*\ , the bispectrum The keyword *bzeroflag* determines whether or not *B0*\ , the bispectrum
components of an atom with no neighbors, are subtracted from components of an atom with no neighbors, are subtracted from
@ -219,13 +228,72 @@ normally only affects compute *sna/atom*\ . However, when
*quadraticflag* is on, it also affects *snad/atom* and *snav/atom*\ . *quadraticflag* is on, it also affects *snad/atom* and *snav/atom*\ .
The keyword *quadraticflag* determines whether or not the The keyword *quadraticflag* determines whether or not the
quadratic analogs to the bispectrum quantities are generated. quadratic combinations of bispectrum quantities are generated.
These are formed by taking the outer product of the vector These are formed by taking the outer product of the vector
of bispectrum components with itself. of bispectrum components with itself.
See section below on output for a See section below on output for a
detailed explanation of the number of quadratic terms and the detailed explanation of the number of quadratic terms and the
ordered in which they are listed. ordered in which they are listed.
The keyword *chem* activates the explicit multi-element variant
of the SNAP bispectrum components. The argument *nelements*
specifies the number of SNAP elements that will be handled.
This is followed by *elementlist*, a list of integers of
length *ntypes*, with values in the range [0, *nelements* ),
which maps each LAMMPS type to one of the SNAP elements.
Note that multiple LAMMPS types can be mapped to the same element,
and some elements may be mapped by no LAMMPS type. However, in typical
use cases (training SNAP potentials) the mapping from LAMMPS types
to elements is one-to-one.
The explicit multi-element variant invoked by the *chem* keyword
partitions the density of neighbors into partial densities
for each chemical element. This is described in detail in the
paper by :ref:`Cusentino et al. <Cusentino2020>`
The bispectrum components are indexed on
ordered triplets of elements:
.. math::
B_{j_1,j_2,j}^{\kappa\lambda\mu} =
\sum_{m_1,m'_1=-j_1}^{j_1}\sum_{m_2,m'_2=-j_2}^{j_2}\sum_{m,m'=-j}^{j} (u^{\mu}_{j,m,m'})^*
H {\scriptscriptstyle \begin{array}{l} {j} {m} {m'} \\
{j_1} {m_1} {m'_1} \\
{j_2} {m_2} {m'_2} \end{array}}
u^{\kappa}_{j_1,m_1,m'_1} u^{\lambda}_{j_2,m_2,m'_2}
where :math:`u^{\mu}_{j,m,m'}` is an expansion coefficient for the partial density of neighbors
of element :math:`\mu`
.. math::
u^{\mu}_{j,m,m'} = w^{self}_{\mu_{i}\mu} U^{j,m,m'}(0,0,0) + \sum_{r_{ii'} < R_{ii'}}{\delta_{\mu\mu_{i'}}f_c(r_{ii'}) w_{\mu_{i'}} U^{j,m,m'}(\theta_0,\theta,\phi)}
where :math:`w^{self}_{\mu_{i}\mu}` is the self-conribution, which is either 1 or 0
(see keyword *wselfallflag* below), :math:`\delta_{\mu\mu_{i'}}` indicates
that the sum is only over neighbor atoms of element :math:`\mu`,
and all other quantities are the same as those appearing in the
original equation for :math:`u^j_{m,m'}` given above.
The keyword *wselfallflag* defines the rule used for the self-contribution.
If *wselfallflag* is on, then :math:`w^{self}_{\mu_{i}\mu}` = 1. If it is
off then :math:`w^{self}_{\mu_{i}\mu}` = 0, except in the case
of :math:`{\mu_{i}=\mu}`, when :math:`w^{self}_{\mu_{i}\mu}` = 1.
When the *chem* keyword is not used, this keyword has no effect.
The keyword *bnormflag* determines whether or not the bispectrum
component :math:`B_{j_1,j_2,j}` is divided by a factor of :math:`2j+1`.
This normalization simplifies force calculations because of the
following symmetry relation
.. math::
\frac{B_{j_1,j_2,j}}{2j+1} = \frac{B_{j,j_2,j_1}}{2j_1+1} = \frac{B_{j_1,j,j_2}}{2j_2+1}
This option is typically used in conjunction with the *chem* keyword,
and LAMMPS will generate a warning if both *chem* and *bnormflag*
are not both set or not both unset.
.. note:: .. note::
If you have a bonded system, then the settings of If you have a bonded system, then the settings of
@ -257,6 +325,8 @@ described by the following piece of python code:
for j in range(j1-j2,min(twojmax,j1+j2)+1,2): for j in range(j1-j2,min(twojmax,j1+j2)+1,2):
if (j>=j1): print j1/2.,j2/2.,j/2. if (j>=j1): print j1/2.,j2/2.,j/2.
For even twojmax = 2(*m*\ -1), :math:`K = m(m+1)(2m+1)/6`, the *m*\ -th pyramidal number. For odd twojmax = 2 *m*\ -1, :math:`K = m(m+1)(m+2)/3`, twice the *m*\ -th tetrahedral number.
.. note:: .. note::
the *diagonal* keyword allowing other possible choices the *diagonal* keyword allowing other possible choices
@ -267,16 +337,15 @@ described by the following piece of python code:
Compute *snad/atom* evaluates a per-atom array. The columns are Compute *snad/atom* evaluates a per-atom array. The columns are
arranged into *ntypes* blocks, listed in order of atom type *I*\ . Each arranged into *ntypes* blocks, listed in order of atom type *I*\ . Each
block contains three sub-blocks corresponding to the *x*\ , *y*\ , and *z* block contains three sub-blocks corresponding to the *x*\ , *y*\ , and *z*
components of the atom position. Each of these sub-blocks contains components of the atom position. Each of these sub-blocks contains *K*
one column for each bispectrum component, the same as for compute columns for the *K* bispectrum components, the same as for compute *sna/atom*
*sna/atom*
Compute *snav/atom* evaluates a per-atom array. The columns are Compute *snav/atom* evaluates a per-atom array. The columns are
arranged into *ntypes* blocks, listed in order of atom type *I*\ . Each arranged into *ntypes* blocks, listed in order of atom type *I*\ . Each
block contains six sub-blocks corresponding to the *xx*\ , *yy*\ , *zz*\ , block contains six sub-blocks corresponding to the *xx*\ , *yy*\ , *zz*\ ,
*yz*\ , *xz*\ , and *xy* components of the virial tensor in Voigt *yz*\ , *xz*\ , and *xy* components of the virial tensor in Voigt
notation. Each of these sub-blocks contains one column for each notation. Each of these sub-blocks contains *K*
bispectrum component, the same as for compute *sna/atom* columns for the *K* bispectrum components, the same as for compute *sna/atom*
Compute *snap* evaluates a global array. Compute *snap* evaluates a global array.
The columns are arranged into The columns are arranged into
@ -312,6 +381,14 @@ of linear terms i.e. linear and quadratic terms are contiguous.
So the nesting order from inside to outside is bispectrum component, So the nesting order from inside to outside is bispectrum component,
linear then quadratic, vector/tensor component, type. linear then quadratic, vector/tensor component, type.
If the *chem* keyword is used, then the data is arranged into :math:`N_{elem}^3`
sub-blocks, each sub-block corresponding to a particular chemical labelling
:math:`\kappa\lambda\mu` with the last label changing fastest.
Each sub-block contains *K* bispectrum components. For the purposes
of handling contributions to force, virial, and quadratic combinations,
these :math:`N_{elem}^3` sub-blocks are treated as a single block
of :math:`K N_{elem}^3` columns.
These values can be accessed by any command that uses per-atom values These values can be accessed by any command that uses per-atom values
from a compute as input. See the :doc:`Howto output <Howto_output>` doc from a compute as input. See the :doc:`Howto output <Howto_output>` doc
page for an overview of LAMMPS output options. page for an overview of LAMMPS output options.
@ -320,7 +397,8 @@ Restrictions
"""""""""""" """"""""""""
These computes are part of the SNAP package. They are only enabled if These computes are part of the SNAP package. They are only enabled if
LAMMPS was built with that package. See the :doc:`Build package <Build_package>` doc page for more info. LAMMPS was built with that package. See the :doc:`Build package <Build_package>`
doc page for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""
@ -332,6 +410,7 @@ Default
The optional keyword defaults are *rmin0* = 0, The optional keyword defaults are *rmin0* = 0,
*switchflag* = 1, *bzeroflag* = 1, *quadraticflag* = 0, *switchflag* = 1, *bzeroflag* = 1, *quadraticflag* = 0,
*bnormflag* = 0, *wselfallflag* = 0
---------- ----------
@ -352,3 +431,12 @@ available at `arXiv:1409.3880 <http://arxiv.org/abs/1409.3880>`_
**(Varshalovich)** Varshalovich, Moskalev, Khersonskii, Quantum Theory **(Varshalovich)** Varshalovich, Moskalev, Khersonskii, Quantum Theory
of Angular Momentum, World Scientific, Singapore (1987). of Angular Momentum, World Scientific, Singapore (1987).
.. _Varshalovich1987:
.. _Mason2009:
**(Mason)** J. K. Mason, Acta Cryst A65, 259 (2009).
.. _Cusentino2020:
**(Cusentino)** Cusentino, Wood, and Thompson, J Phys Chem A, xxx, xxxxx, (2020)

76
doc/src/pair_snap.rst
View File

@ -24,27 +24,30 @@ Examples
Description Description
""""""""""" """""""""""
Pair style *snap* computes interactions using the spectral Pair style *snap* defines the spectral
neighbor analysis potential (SNAP) :ref:`(Thompson) <Thompson20142>`. neighbor analysis potential (SNAP), a machine-learning
interatomic potential :ref:`(Thompson) <Thompson20142>`.
Like the GAP framework of Bartok et al. :ref:`(Bartok2010) <Bartok20102>`, Like the GAP framework of Bartok et al. :ref:`(Bartok2010) <Bartok20102>`,
:ref:`(Bartok2013) <Bartok2013>` which uses bispectrum components SNAP uses bispectrum components
to characterize the local neighborhood of each atom to characterize the local neighborhood of each atom
in a very general way. The mathematical definition of the in a very general way. The mathematical definition of the
bispectrum calculation used by SNAP is identical bispectrum calculation and its derivatives w.r.t. atom positions
to that used by :doc:`compute sna/atom <compute_sna_atom>`. is identical to that used by :doc:`compute snap <compute_sna_atom>`,
which is used to fit SNAP potentials to *ab initio* energy, force,
and stress data.
In SNAP, the total energy is decomposed into a sum over In SNAP, the total energy is decomposed into a sum over
atom energies. The energy of atom *i* is atom energies. The energy of atom *i* is
expressed as a weighted sum over bispectrum components. expressed as a weighted sum over bispectrum components.
.. math:: .. math::
E^i_{SNAP}(B_1^i,...,B_K^i) = \beta^{\alpha_i}_0 + \sum_{k=1}^K \beta_k^{\alpha_i} B_k^i E^i_{SNAP}(B_1^i,...,B_K^i) = \beta^{\mu_i}_0 + \sum_{k=1}^K \beta_k^{\mu_i} B_k^i
where :math:`B_k^i` is the *k*\ -th bispectrum component of atom *i*\ , where :math:`B_k^i` is the *k*\ -th bispectrum component of atom *i*\ ,
and :math:`\beta_k^{\alpha_i}` is the corresponding linear coefficient and :math:`\beta_k^{\mu_i}` is the corresponding linear coefficient
that depends on :math:\alpha_i`, the SNAP element of atom *i*\ . The that depends on :math:`\mu_i`, the SNAP element of atom *i*\ . The
number of bispectrum components used and their definitions number of bispectrum components used and their definitions
depend on the value of *twojmax* depend on the value of *twojmax* and other parameters
defined in the SNAP parameter file described below. defined in the SNAP parameter file described below.
The bispectrum calculation is described in more detail The bispectrum calculation is described in more detail
in :doc:`compute sna/atom <compute_sna_atom>`. in :doc:`compute sna/atom <compute_sna_atom>`.
@ -136,17 +139,51 @@ The SNAP parameter file can contain blank and comment lines (start
with #) anywhere. Each non-blank non-comment line must contain one with #) anywhere. Each non-blank non-comment line must contain one
keyword/value pair. The required keywords are *rcutfac* and keyword/value pair. The required keywords are *rcutfac* and
*twojmax*\ . Optional keywords are *rfac0*\ , *rmin0*\ , *twojmax*\ . Optional keywords are *rfac0*\ , *rmin0*\ ,
*switchflag*\ , *bzeroflag*\, and *chunksize*\. *switchflag*\ , *bzeroflag*\ , *quadraticflag*\ , *chemflag*\ ,
*bnormflag*\ , *wselfallflag*\ , and *chunksize*\ .
The default values for these keywords are The default values for these keywords are
* *rfac0* = 0.99363 * *rfac0* = 0.99363
* *rmin0* = 0.0 * *rmin0* = 0.0
* *switchflag* = 0 * *switchflag* = 1
* *bzeroflag* = 1 * *bzeroflag* = 1
* *quadraticflag* = 1 * *quadraticflag* = 0
* *chemflag* = 0
* *bnormflag* = 0
* *wselfallflag* = 0
* *chunksize* = 2000 * *chunksize* = 2000
If *quadraticflag* is set to 1, then the SNAP energy expression includes additional quadratic terms
that have been shown to increase the overall accuracy of the potential without much increase
in computational cost :ref:`(Wood) <Wood20182>`.
.. math::
E^i_{SNAP}(\mathbf{B}^i) = \beta^{\mu_i}_0 + \boldsymbol{\beta}^{\mu_i} \cdot \mathbf{B}_i + \frac{1}{2}\mathbf{B}^t_i \cdot \boldsymbol{\alpha}^{\mu_i} \cdot \mathbf{B}_i
where :math:`\mathbf{B}_i` is the *K*-vector of bispectrum components,
:math:`\boldsymbol{\beta}^{\mu_i}` is the *K*-vector of linear coefficients
for element :math:`\mu_i`, and :math:`\boldsymbol{\alpha}^{\mu_i}`
is the symmetric *K* by *K* matrix of quadratic coefficients.
The SNAP element file should contain *K*\ (\ *K*\ +1)/2 additional coefficients
for each element, the upper-triangular elements of :math:`\boldsymbol{\alpha}^{\mu_i}`.
If *chemflag* is set to 1, then the energy expression is written in terms of explicit multi-element bispectrum
components indexed on ordered triplets of elements, which has been shown to increase the ability of the SNAP
potential to capture energy differences in chemically complex systems,
at the expense of a significant increase in computational cost :ref:`(Cusentino) <Cusentino20202>`.
.. math::
E^i_{SNAP}(\mathbf{B}^i) = \beta^{\mu_i}_0 + \sum_{\kappa,\lambda,\mu} \boldsymbol{\beta}^{\kappa\lambda\mu}_{\mu_i} \cdot \mathbf{B}^{\kappa\lambda\mu}_i
where :math:`\mathbf{B}^{\kappa\lambda\mu}_i` is the *K*-vector of bispectrum components
for neighbors of elements :math:`\kappa`, :math:`\lambda`, and :math:`\mu` and
:math:`\boldsymbol{\beta}^{\kappa\lambda\mu}_{\mu_i}` is the corresponding *K*-vector
of linear coefficients for element :math:`\mu_i`. The SNAP element file should contain
a total of :math:`K N_{elem}^3` coefficients for each of the :math:`N_{elem}` elements.
The keyword *chunksize* is only applicable when using the The keyword *chunksize* is only applicable when using the
pair style *snap* with the KOKKOS package and is ignored otherwise. pair style *snap* with the KOKKOS package and is ignored otherwise.
This keyword controls This keyword controls
@ -159,10 +196,6 @@ into two passes.
Detailed definitions for all the other keywords Detailed definitions for all the other keywords
are given on the :doc:`compute sna/atom <compute_sna_atom>` doc page. are given on the :doc:`compute sna/atom <compute_sna_atom>` doc page.
If *quadraticflag* is set to 1, then the SNAP energy expression includes the quadratic term, 0.5\*B\^t.alpha.B, where alpha is a symmetric *K* by *K* matrix.
The SNAP element file should contain *K*\ (\ *K*\ +1)/2 additional coefficients
for each element, the upper-triangular elements of alpha.
.. note:: .. note::
The previously used *diagonalstyle* keyword was removed in 2019, The previously used *diagonalstyle* keyword was removed in 2019,
@ -221,7 +254,8 @@ Related commands
:doc:`compute sna/atom <compute_sna_atom>`, :doc:`compute sna/atom <compute_sna_atom>`,
:doc:`compute snad/atom <compute_sna_atom>`, :doc:`compute snad/atom <compute_sna_atom>`,
:doc:`compute snav/atom <compute_sna_atom>` :doc:`compute snav/atom <compute_sna_atom>`,
:doc:`compute snap <compute_sna_atom>`
**Default:** none **Default:** none
@ -235,6 +269,10 @@ Related commands
**(Bartok2010)** Bartok, Payne, Risi, Csanyi, Phys Rev Lett, 104, 136403 (2010). **(Bartok2010)** Bartok, Payne, Risi, Csanyi, Phys Rev Lett, 104, 136403 (2010).
.. _Bartok2013: .. _Wood20182:
**(Bartok2013)** Bartok, Gillan, Manby, Csanyi, Phys Rev B 87, 184115 (2013). **(Wood)** Wood and Thompson, J Chem Phys, 148, 241721, (2018)
.. _Cusentino20202:
**(Cusentino)** Cusentino, Wood, and Thompson, J Phys Chem A, xxx, xxxxx, (2020)

21
examples/snap/InP_Cusentino_PRB2020.snap
View File

@ -1,21 +0,0 @@
# DATE: 2019-09-18 CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Wood, M.A. Cusentino, B.D. Wirth, and A.P. Thompson, "Data-driven material models for atomistic simulation", Physical Review B 99, 184305 (2019)
# Definition of SNAP+ZBL potential.
set type 1 charge 1e-08
set type 2 charge -1e-08
variable zblcutinner equal 4
variable zblcutouter equal 4.2
variable zblz1 equal 49
variable zblz2 equal 15
variable rcoul equal 10
# Specify hybrid with SNAP, ZBL, and long-range Coulomb
pair_style hybrid/overlay coul/long ${rcoul} &
zbl ${zblcutinner} ${zblcutouter} &
snap
pair_coeff * * coul/long
pair_coeff 1 1 zbl ${zblz1} ${zblz1}
pair_coeff 1 2 zbl ${zblz1} ${zblz2}
pair_coeff 2 2 zbl ${zblz2} ${zblz2}
pair_coeff * * snap InP_Cusentino_PRB2020.snapcoeff InP_Cusentino_PRB2020.snapparam In P
kspace_style ewald 1.0e-5

13
examples/snap/InP_Cusentino_PRB2020.snapparam
View File

@ -1,13 +0,0 @@
# required
rcutfac 1.0
twojmax 6
# optional
rfac0 0.99363
rmin0 0.0
bzeroflag 1
quadraticflag 0
wselfallflag 1
alloyflag 1

1
examples/snap/InP_JCPA2020.snap Symbolic link
View File

@ -0,0 +1 @@
../../potentials/InP_JCPA2020.snap

1
examples/snap/InP_JCPA2020.snapcoeff Symbolic link
View File

@ -0,0 +1 @@
../../potentials/InP_JCPA2020.snapcoeff

1
examples/snap/InP_JCPA2020.snapparam Symbolic link
View File

@ -0,0 +1 @@
../../potentials/InP_JCPA2020.snapparam

3
examples/snap/in.snap.InP.PRB2020 → examples/snap/in.snap.InP.JCPA2020
View File

@ -6,7 +6,6 @@ variable nsteps index 100
variable nrep equal 4 variable nrep equal 4
variable a equal 5.83 variable a equal 5.83
units metal units metal
atom_style charge
# generate the box and atom positions using a FCC lattice # generate the box and atom positions using a FCC lattice
@ -26,7 +25,7 @@ mass 2 30.98
# choose potential # choose potential
include InP_Cusentino_PRB2020.snap include InP_JCPA2020.snap
# Setup output # Setup output

156
examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.1 Normal file
View File

@ -0,0 +1,156 @@
LAMMPS (2 Jun 2020)
# Demonstrate SNAP InP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 5.83
units metal
# generate the box and atom positions using a FCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice diamond $a
lattice diamond 5.83
Lattice spacing in x,y,z = 5.83 5.83 5.83
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 2 box
Created orthogonal box = (0.0 0.0 0.0) to (23.32 23.32 23.32)
1 by 1 by 1 MPI processor grid
create_atoms 1 box basis 5 2 basis 6 2 basis 7 2 basis 8 2
Created 512 atoms
create_atoms CPU = 0.000 seconds
mass 1 114.76
mass 2 30.98
# choose potential
include InP_JCPA2020.snap
# DATE: 2020-06-01 CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020)
# Definition of SNAP+ZBL potential.
variable zblcutinner equal 4
variable zblcutouter equal 4.2
variable zblz1 equal 49
variable zblz2 equal 15
# Specify hybrid with SNAP and ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} snap
pair_style hybrid/overlay zbl 4 ${zblcutouter} snap
pair_style hybrid/overlay zbl 4 4.2 snap
pair_coeff 1 1 zbl ${zblz1} ${zblz1}
pair_coeff 1 1 zbl 49 ${zblz1}
pair_coeff 1 1 zbl 49 49
pair_coeff 1 2 zbl ${zblz1} ${zblz2}
pair_coeff 1 2 zbl 49 ${zblz2}
pair_coeff 1 2 zbl 49 15
pair_coeff 2 2 zbl ${zblz2} ${zblz2}
pair_coeff 2 2 zbl 15 ${zblz2}
pair_coeff 2 2 zbl 15 15
pair_coeff * * snap InP_JCPA2020.snapcoeff InP_JCPA2020.snapparam In P
Reading potential file InP_JCPA2020.snapcoeff with DATE: 2020-06-01
SNAP Element = In, Radius 3.81205, Weight 1
SNAP Element = P, Radius 3.82945, Weight 0.929316
Reading potential file InP_JCPA2020.snapparam with DATE: 2020-06-01
SNAP keyword rcutfac 1.0
SNAP keyword twojmax 6
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0.0
SNAP keyword bzeroflag 1
SNAP keyword quadraticflag 0
SNAP keyword wselfallflag 1
SNAP keyword chemflag 1
SNAP keyword bnormflag 1
# Setup output
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8.6589
ghost atom cutoff = 8.6589
binsize = 4.32945, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair snap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 6.027 | 6.027 | 6.027 Mbytes
Step Temp E_pair E_mol TotEng Press
0 300 -3.4805794 0 -3.4418771 1353.5968
10 286.42264 -3.4788274 0 -3.4418767 1586.4881
20 250.14148 -3.4741459 0 -3.4418757 2219.0344
30 202.52417 -3.4680017 0 -3.4418745 2982.7272
40 157.39113 -3.4621782 0 -3.4418735 3631.0936
50 126.7004 -3.4582183 0 -3.441873 4053.7725
60 117.00722 -3.4569679 0 -3.441873 4204.9542
70 127.65968 -3.4583427 0 -3.4418736 4106.3112
80 151.50217 -3.4614195 0 -3.4418745 3840.7205
90 177.67607 -3.464797 0 -3.4418754 3527.9794
100 195.30359 -3.4670717 0 -3.4418761 3300.3795
Loop time of 18.0803 on 1 procs for 100 steps with 512 atoms
Performance: 0.239 ns/day, 100.446 hours/ns, 5.531 timesteps/s
99.8% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 18.078 | 18.078 | 18.078 | 0.0 | 99.99
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.000979 | 0.000979 | 0.000979 | 0.0 | 0.01
Output | 0.000243 | 0.000243 | 0.000243 | 0.0 | 0.00
Modify | 0.000591 | 0.000591 | 0.000591 | 0.0 | 0.00
Other | | 0.000469 | | | 0.00
Nlocal: 512 ave 512 max 512 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 1959 ave 1959 max 1959 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 31232 ave 31232 max 31232 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 62464 ave 62464 max 62464 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 62464
Ave neighs/atom = 122
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:18

156
examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.4 Normal file
View File

@ -0,0 +1,156 @@
LAMMPS (2 Jun 2020)
# Demonstrate SNAP InP potential
# Initialize simulation
variable nsteps index 100
variable nrep equal 4
variable a equal 5.83
units metal
# generate the box and atom positions using a FCC lattice
variable nx equal ${nrep}
variable nx equal 4
variable ny equal ${nrep}
variable ny equal 4
variable nz equal ${nrep}
variable nz equal 4
boundary p p p
lattice diamond $a
lattice diamond 5.83
Lattice spacing in x,y,z = 5.83 5.83 5.83
region box block 0 ${nx} 0 ${ny} 0 ${nz}
region box block 0 4 0 ${ny} 0 ${nz}
region box block 0 4 0 4 0 ${nz}
region box block 0 4 0 4 0 4
create_box 2 box
Created orthogonal box = (0.0 0.0 0.0) to (23.32 23.32 23.32)
1 by 2 by 2 MPI processor grid
create_atoms 1 box basis 5 2 basis 6 2 basis 7 2 basis 8 2
Created 512 atoms
create_atoms CPU = 0.001 seconds
mass 1 114.76
mass 2 30.98
# choose potential
include InP_JCPA2020.snap
# DATE: 2020-06-01 CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020)
# Definition of SNAP+ZBL potential.
variable zblcutinner equal 4
variable zblcutouter equal 4.2
variable zblz1 equal 49
variable zblz2 equal 15
# Specify hybrid with SNAP and ZBL
pair_style hybrid/overlay zbl ${zblcutinner} ${zblcutouter} snap
pair_style hybrid/overlay zbl 4 ${zblcutouter} snap
pair_style hybrid/overlay zbl 4 4.2 snap
pair_coeff 1 1 zbl ${zblz1} ${zblz1}
pair_coeff 1 1 zbl 49 ${zblz1}
pair_coeff 1 1 zbl 49 49
pair_coeff 1 2 zbl ${zblz1} ${zblz2}
pair_coeff 1 2 zbl 49 ${zblz2}
pair_coeff 1 2 zbl 49 15
pair_coeff 2 2 zbl ${zblz2} ${zblz2}
pair_coeff 2 2 zbl 15 ${zblz2}
pair_coeff 2 2 zbl 15 15
pair_coeff * * snap InP_JCPA2020.snapcoeff InP_JCPA2020.snapparam In P
Reading potential file InP_JCPA2020.snapcoeff with DATE: 2020-06-01
SNAP Element = In, Radius 3.81205, Weight 1
SNAP Element = P, Radius 3.82945, Weight 0.929316
Reading potential file InP_JCPA2020.snapparam with DATE: 2020-06-01
SNAP keyword rcutfac 1.0
SNAP keyword twojmax 6
SNAP keyword rfac0 0.99363
SNAP keyword rmin0 0.0
SNAP keyword bzeroflag 1
SNAP keyword quadraticflag 0
SNAP keyword wselfallflag 1
SNAP keyword chemflag 1
SNAP keyword bnormflag 1
# Setup output
thermo 10
thermo_modify norm yes
# Set up NVE run
timestep 0.5e-3
neighbor 1.0 bin
neigh_modify once no every 1 delay 0 check yes
# Run MD
velocity all create 300.0 4928459
fix 1 all nve
run ${nsteps}
run 100
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 8.6589
ghost atom cutoff = 8.6589
binsize = 4.32945, bins = 6 6 6
2 neighbor lists, perpetual/occasional/extra = 2 0 0
(1) pair zbl, perpetual, half/full from (2)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair snap, perpetual
attributes: full, newton on
pair build: full/bin/atomonly
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.587 | 4.587 | 4.587 Mbytes
Step Temp E_pair E_mol TotEng Press
0 300 -3.4805794 0 -3.4418771 1353.5968
10 286.58246 -3.478848 0 -3.4418767 1582.995
20 250.70996 -3.4742192 0 -3.4418757 2207.7507
30 203.58199 -3.4681382 0 -3.4418746 2968.5206
40 158.84622 -3.462366 0 -3.4418736 3619.0285
50 128.30488 -3.4584254 0 -3.4418731 4047.173
60 118.40349 -3.4571481 0 -3.4418731 4203.3421
70 128.48973 -3.4584499 0 -3.4418737 4109.0296
80 151.54241 -3.4614247 0 -3.4418746 3847.4617
90 176.92084 -3.4646996 0 -3.4418755 3548.7811
100 193.9555 -3.4668978 0 -3.4418761 3342.8083
Loop time of 4.99339 on 4 procs for 100 steps with 512 atoms
Performance: 0.865 ns/day, 27.741 hours/ns, 20.026 timesteps/s
99.5% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 4.8898 | 4.907 | 4.9329 | 0.8 | 98.27
Neigh | 0 | 0 | 0 | 0.0 | 0.00
Comm | 0.058815 | 0.084739 | 0.1019 | 6.1 | 1.70
Output | 0.000252 | 0.00038775 | 0.000777 | 0.0 | 0.01
Modify | 0.000262 | 0.00026675 | 0.000272 | 0.0 | 0.01
Other | | 0.001039 | | | 0.02
Nlocal: 128 ave 128 max 128 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Nghost: 1099 ave 1099 max 1099 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Neighs: 7808 ave 7808 max 7808 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 15616 ave 15616 max 15616 min
Histogram: 4 0 0 0 0 0 0 0 0 0
Total # of neighbors = 62464
Ave neighs/atom = 122
Neighbor list builds = 0
Dangerous builds = 0
Total wall time: 0:00:05

18
potentials/InP_JCPA2020.snap Normal file
View File

@ -0,0 +1,18 @@
# DATE: 2020-06-01 CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020)
# Definition of SNAP+ZBL potential.
variable zblcutinner equal 4
variable zblcutouter equal 4.2
variable zblz1 equal 49
variable zblz2 equal 15
# Specify hybrid with SNAP and ZBL
pair_style hybrid/overlay &
zbl ${zblcutinner} ${zblcutouter} &
snap
pair_coeff 1 1 zbl ${zblz1} ${zblz1}
pair_coeff 1 2 zbl ${zblz1} ${zblz2}
pair_coeff 2 2 zbl ${zblz2} ${zblz2}
pair_coeff * * snap InP_JCPA2020.snapcoeff InP_JCPA2020.snapparam In P

2
examples/snap/InP_Cusentino_PRB2020.snapcoeff → potentials/InP_JCPA2020.snapcoeff
View File

@ -1,4 +1,4 @@
# LAMMPS SNAP coefficients for InP # DATE: 2020-06-01 CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020)
2 241 2 241
In 3.81205 1 In 3.81205 1

15
potentials/InP_JCPA2020.snapparam Normal file
View File

@ -0,0 +1,15 @@
# DATE: 2020-06-01 CONTRIBUTOR: Mary Alice Cusentino mcusent@sandia.gov CITATION: M.A. Cusentino, M. A. Wood, and A.P. Thompson, "Explicit Multi-element Extension of the Spectral Neighbor Analysis Potential for Chemically Complex Systems", J. Phys. Chem. A, xxxxxx (2020)
# required
rcutfac 1.0
twojmax 6
# optional
rfac0 0.99363
rmin0 0.0
bzeroflag 1
quadraticflag 0
wselfallflag 1
chemflag 1
bnormflag 1

4
src/KOKKOS/pair_snap_kokkos_impl.h
View File

@ -515,7 +515,7 @@ void PairSNAPKokkos<DeviceType>::coeff(int narg, char **arg)
Kokkos::deep_copy(d_map,h_map); Kokkos::deep_copy(d_map,h_map);
snaKK = SNAKokkos<DeviceType>(rfac0,twojmax, snaKK = SNAKokkos<DeviceType>(rfac0,twojmax,
rmin0,switchflag,bzeroflag,alloyflag,wselfallflag,nelements); rmin0,switchflag,bzeroflag,chemflag,bnormflag,wselfallflag,nelements);
snaKK.grow_rij(0,0); snaKK.grow_rij(0,0);
snaKK.init(); snaKK.init();
} }
@ -584,7 +584,7 @@ void PairSNAPKokkos<DeviceType>::operator() (TagPairSNAPComputeNeigh,const typen
my_sna.inside(ii,offset) = j; my_sna.inside(ii,offset) = j;
my_sna.wj(ii,offset) = d_wjelem[elem_j]; my_sna.wj(ii,offset) = d_wjelem[elem_j];
my_sna.rcutij(ii,offset) = (radi + d_radelem[elem_j])*rcutfac; my_sna.rcutij(ii,offset) = (radi + d_radelem[elem_j])*rcutfac;
if (alloyflag) if (chemflag)
my_sna.element(ii,offset) = elem_j; my_sna.element(ii,offset) = elem_j;
else else
my_sna.element(ii,offset) = 0; my_sna.element(ii,offset) = 0;

2
src/KOKKOS/sna_kokkos.h
View File

@ -210,7 +210,7 @@ inline
int switch_flag; int switch_flag;
// Chem snap flags // Chem snap flags
int alloy_flag; int chem_flag;
int bnorm_flag; int bnorm_flag;
int nelements; int nelements;
int ndoubles; int ndoubles;

14
src/KOKKOS/sna_kokkos_impl.h
View File

@ -29,7 +29,7 @@ template<class DeviceType>
inline inline
SNAKokkos<DeviceType>::SNAKokkos(double rfac0_in, SNAKokkos<DeviceType>::SNAKokkos(double rfac0_in,
int twojmax_in, double rmin0_in, int switch_flag_in, int bzero_flag_in, int twojmax_in, double rmin0_in, int switch_flag_in, int bzero_flag_in,
int alloy_flag_in, int wselfall_flag_in, int nelements_in) int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in)
{ {
wself = 1.0; wself = 1.0;
@ -38,13 +38,13 @@ SNAKokkos<DeviceType>::SNAKokkos(double rfac0_in,
switch_flag = switch_flag_in; switch_flag = switch_flag_in;
bzero_flag = bzero_flag_in; bzero_flag = bzero_flag_in;
alloy_flag = alloy_flag_in; chem_flag = chem_flag_in;
wselfall_flag = wselfall_flag_in; if (chem_flag)
if (alloy_flag)
nelements = nelements_in; nelements = nelements_in;
else else
nelements = 1; nelements = 1;
bnorm_flag = alloy_flag_in; bnorm_flag = bnorm_flag_in;
wselfall_flag = wselfall_flag_in;
twojmax = twojmax_in; twojmax = twojmax_in;
@ -283,7 +283,7 @@ void SNAKokkos<DeviceType>::pre_ui(const typename Kokkos::TeamPolicy<DeviceType>
// if m is on the "diagonal", initialize it with the self energy. // if m is on the "diagonal", initialize it with the self energy.
// Otherwise zero it out // Otherwise zero it out
SNAcomplex init = {0., 0.}; SNAcomplex init = {0., 0.};
if (m % (j+2) == 0 && (!alloy_flag || ielem == jelem || wselfall_flag)) { init = {wself, 0.0}; } //need to map iatom to element if (m % (j+2) == 0 && (!chem_flag || ielem == jelem || wselfall_flag)) { init = {wself, 0.0}; } //need to map iatom to element
ulisttot(jelem*idxu_max+jjup, iatom) = init; ulisttot(jelem*idxu_max+jjup, iatom) = init;
}); });
@ -1620,7 +1620,7 @@ int SNAKokkos<DeviceType>::compute_ncoeff()
ndoubles = nelements*nelements; ndoubles = nelements*nelements;
ntriples = nelements*nelements*nelements; ntriples = nelements*nelements*nelements;
if (alloy_flag) ncount *= ntriples; if (chem_flag) ncount *= ntriples;
return ncount; return ncount;
} }

24
src/SNAP/compute_sna_atom.cpp
View File

@ -34,7 +34,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
radelem(NULL), wjelem(NULL) radelem(NULL), wjelem(NULL)
{ {
double rmin0, rfac0; double rmin0, rfac0;
int twojmax, switchflag, bzeroflag, bnormflag; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
radelem = NULL; radelem = NULL;
wjelem = NULL; wjelem = NULL;
@ -50,7 +50,8 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
bzeroflag = 1; bzeroflag = 1;
bnormflag = 0; bnormflag = 0;
quadraticflag = 0; quadraticflag = 0;
alloyflag = 0; chemflag = 0;
bnormflag = 0;
wselfallflag = 0; wselfallflag = 0;
nelements = 1; nelements = 1;
@ -108,22 +109,25 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Illegal compute sna/atom command"); error->all(FLERR,"Illegal compute sna/atom command");
quadraticflag = atoi(arg[iarg+1]); quadraticflag = atoi(arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"alloy") == 0) { } else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg) if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute sna/atom command"); error->all(FLERR,"Illegal compute sna/atom command");
alloyflag = 1; chemflag = 1;
bnormflag = alloyflag;
memory->create(map,ntypes+1,"compute_sna_atom:map"); memory->create(map,ntypes+1,"compute_sna_atom:map");
nelements = force->inumeric(FLERR,arg[iarg+1]); nelements = force->inumeric(FLERR,arg[iarg+1]);
for(int i = 0; i < ntypes; i++) { for(int i = 0; i < ntypes; i++) {
int jelem = force->inumeric(FLERR,arg[iarg+2+i]); int jelem = force->inumeric(FLERR,arg[iarg+2+i]);
printf("%d %d %d %d\n",ntypes,nelements,i,jelem);
if (jelem < 0 || jelem >= nelements) if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute sna/atom command"); error->all(FLERR,"Illegal compute sna/atom command");
map[i+1] = jelem; map[i+1] = jelem;
} }
iarg += 2+ntypes; iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"wselfall") == 0) { } else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute sna/atom command"); error->all(FLERR,"Illegal compute sna/atom command");
wselfallflag = atoi(arg[iarg+1]); wselfallflag = atoi(arg[iarg+1]);
@ -133,7 +137,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) :
snaptr = new SNA(lmp, rfac0, twojmax, snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag, rmin0, switchflag, bzeroflag,
alloyflag, wselfallflag, nelements); chemflag, bnormflag, wselfallflag, nelements);
ncoeff = snaptr->ncoeff; ncoeff = snaptr->ncoeff;
size_peratom_cols = ncoeff; size_peratom_cols = ncoeff;
@ -229,7 +233,7 @@ void ComputeSNAAtom::compute_peratom()
const double ztmp = x[i][2]; const double ztmp = x[i][2];
const int itype = type[i]; const int itype = type[i];
int ielem = 0; int ielem = 0;
if (alloyflag) if (chemflag)
ielem = map[itype]; ielem = map[itype];
const double radi = radelem[itype]; const double radi = radelem[itype];
const int* const jlist = firstneigh[i]; const int* const jlist = firstneigh[i];
@ -254,7 +258,7 @@ void ComputeSNAAtom::compute_peratom()
const double rsq = delx*delx + dely*dely + delz*delz; const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j]; int jtype = type[j];
int jelem = 0; int jelem = 0;
if (alloyflag) if (chemflag)
int jelem = map[jtype]; int jelem = map[jtype];
if (rsq < cutsq[itype][jtype] && rsq>1e-20) { if (rsq < cutsq[itype][jtype] && rsq>1e-20) {
snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][0] = delx;

2
src/SNAP/compute_sna_atom.h
View File

@ -43,7 +43,7 @@ class ComputeSNAAtom : public Compute {
double *radelem; double *radelem;
double *wjelem; double *wjelem;
int * map; // map types to [0,nelements) int * map; // map types to [0,nelements)
int nelements, alloyflag, wselfallflag; int nelements, chemflag;
class SNA* snaptr; class SNA* snaptr;
double cutmax; double cutmax;
int quadraticflag; int quadraticflag;

24
src/SNAP/compute_snad_atom.cpp
View File

@ -34,7 +34,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
radelem(NULL), wjelem(NULL) radelem(NULL), wjelem(NULL)
{ {
double rfac0, rmin0; double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
radelem = NULL; radelem = NULL;
wjelem = NULL; wjelem = NULL;
@ -50,7 +50,8 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
bzeroflag = 1; bzeroflag = 1;
bnormflag = 0; bnormflag = 0;
quadraticflag = 0; quadraticflag = 0;
alloyflag = 0; chemflag = 0;
bnormflag = 0;
wselfallflag = 0; wselfallflag = 0;
nelements = 1; nelements = 1;
@ -106,22 +107,25 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Illegal compute snad/atom command"); error->all(FLERR,"Illegal compute snad/atom command");
quadraticflag = atoi(arg[iarg+1]); quadraticflag = atoi(arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"alloy") == 0) { } else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg) if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute snad/atom command"); error->all(FLERR,"Illegal compute snad/atom command");
alloyflag = 1; chemflag = 1;
bnormflag = alloyflag;
memory->create(map,ntypes+1,"compute_snad_atom:map"); memory->create(map,ntypes+1,"compute_snad_atom:map");
nelements = force->inumeric(FLERR,arg[iarg+1]); nelements = force->inumeric(FLERR,arg[iarg+1]);
for(int i = 0; i < ntypes; i++) { for(int i = 0; i < ntypes; i++) {
int jelem = force->inumeric(FLERR,arg[iarg+2+i]); int jelem = force->inumeric(FLERR,arg[iarg+2+i]);
printf("%d %d %d %d\n",ntypes,nelements,i,jelem);
if (jelem < 0 || jelem >= nelements) if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snad/atom command"); error->all(FLERR,"Illegal compute snad/atom command");
map[i+1] = jelem; map[i+1] = jelem;
} }
iarg += 2+ntypes; iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"wselfall") == 0) { } else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snad/atom command"); error->all(FLERR,"Illegal compute snad/atom command");
wselfallflag = atoi(arg[iarg+1]); wselfallflag = atoi(arg[iarg+1]);
@ -131,7 +135,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) :
snaptr = new SNA(lmp, rfac0, twojmax, snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag, rmin0, switchflag, bzeroflag,
alloyflag, wselfallflag, nelements); chemflag, bnormflag, wselfallflag, nelements);
ncoeff = snaptr->ncoeff; ncoeff = snaptr->ncoeff;
nperdim = ncoeff; nperdim = ncoeff;
@ -241,7 +245,7 @@ void ComputeSNADAtom::compute_peratom()
const double ztmp = x[i][2]; const double ztmp = x[i][2];
const int itype = type[i]; const int itype = type[i];
int ielem = 0; int ielem = 0;
if (alloyflag) if (chemflag)
ielem = map[itype]; ielem = map[itype];
const double radi = radelem[itype]; const double radi = radelem[itype];
const int* const jlist = firstneigh[i]; const int* const jlist = firstneigh[i];
@ -272,7 +276,7 @@ void ComputeSNADAtom::compute_peratom()
const double rsq = delx*delx + dely*dely + delz*delz; const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j]; int jtype = type[j];
int jelem = 0; int jelem = 0;
if (alloyflag) if (chemflag)
jelem = map[jtype]; jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][0] = delx;

2
src/SNAP/compute_snad_atom.h
View File

@ -45,7 +45,7 @@ class ComputeSNADAtom : public Compute {
double *radelem; double *radelem;
double *wjelem; double *wjelem;
int *map; // map types to [0,nelements) int *map; // map types to [0,nelements)
int nelements, alloyflag, wselfallflag; int nelements, chemflag;
class SNA* snaptr; class SNA* snaptr;
double cutmax; double cutmax;
int quadraticflag; int quadraticflag;

23
src/SNAP/compute_snap.cpp
View File

@ -41,7 +41,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
extarray = 0; extarray = 0;
double rfac0, rmin0; double rfac0, rmin0;
int twojmax, switchflag, bzeroflag; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
radelem = NULL; radelem = NULL;
wjelem = NULL; wjelem = NULL;
@ -56,7 +56,8 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
switchflag = 1; switchflag = 1;
bzeroflag = 1; bzeroflag = 1;
quadraticflag = 0; quadraticflag = 0;
alloyflag = 0; chemflag = 0;
bnormflag = 0;
wselfallflag = 0; wselfallflag = 0;
nelements = 1; nelements = 1;
@ -112,20 +113,24 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Illegal compute snap command"); error->all(FLERR,"Illegal compute snap command");
quadraticflag = atoi(arg[iarg+1]); quadraticflag = atoi(arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"alloyflag") == 0) { } else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command"); error->all(FLERR,"Illegal compute snap command");
alloyflag = 1; chemflag = 1;
memory->create(map,ntypes+1,"compute_snap:map"); memory->create(map,ntypes+1,"compute_snap:map");
nelements = force->inumeric(FLERR,arg[iarg+1]); nelements = force->inumeric(FLERR,arg[iarg+1]);
for(int i = 0; i < ntypes; i++) { for(int i = 0; i < ntypes; i++) {
int jelem = force->inumeric(FLERR,arg[iarg+2+i]); int jelem = force->inumeric(FLERR,arg[iarg+2+i]);
if (screen && comm->me==0) fprintf(screen, "%d %d %d %d\n",ntypes,nelements,i,jelem);
if (jelem < 0 || jelem >= nelements) if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snap command"); error->all(FLERR,"Illegal compute snap command");
map[i+1] = jelem; map[i+1] = jelem;
} }
iarg += 2+ntypes; iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) { } else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snap command"); error->all(FLERR,"Illegal compute snap command");
@ -136,7 +141,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) :
snaptr = new SNA(lmp, rfac0, twojmax, snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag, rmin0, switchflag, bzeroflag,
alloyflag, wselfallflag, nelements); chemflag, bnormflag, wselfallflag, nelements);
ncoeff = snaptr->ncoeff; ncoeff = snaptr->ncoeff;
nperdim = ncoeff; nperdim = ncoeff;
@ -168,7 +173,7 @@ ComputeSnap::~ComputeSnap()
memory->destroy(cutsq); memory->destroy(cutsq);
delete snaptr; delete snaptr;
if (alloyflag) memory->destroy(map); if (chemflag) memory->destroy(map);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -300,7 +305,7 @@ void ComputeSnap::compute_array()
const double ztmp = x[i][2]; const double ztmp = x[i][2];
const int itype = type[i]; const int itype = type[i];
int ielem = 0; int ielem = 0;
if (alloyflag) if (chemflag)
ielem = map[itype]; ielem = map[itype];
const double radi = radelem[itype]; const double radi = radelem[itype];
const int* const jlist = firstneigh[i]; const int* const jlist = firstneigh[i];
@ -328,7 +333,7 @@ void ComputeSnap::compute_array()
const double rsq = delx*delx + dely*dely + delz*delz; const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j]; int jtype = type[j];
int jelem = 0; int jelem = 0;
if (alloyflag) if (chemflag)
jelem = map[jtype]; jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][0] = delx;

2
src/SNAP/compute_snap.h
View File

@ -45,7 +45,7 @@ class ComputeSnap : public Compute {
double *radelem; double *radelem;
double *wjelem; double *wjelem;
int *map; // map types to [0,nelements) int *map; // map types to [0,nelements)
int nelements, alloyflag, wselfallflag; int nelements, chemflag;
class SNA* snaptr; class SNA* snaptr;
double cutmax; double cutmax;
int quadraticflag; int quadraticflag;

24
src/SNAP/compute_snav_atom.cpp
View File

@ -33,7 +33,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
radelem(NULL), wjelem(NULL) radelem(NULL), wjelem(NULL)
{ {
double rfac0, rmin0; double rfac0, rmin0;
int twojmax, switchflag, bzeroflag, bnormflag; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag;
radelem = NULL; radelem = NULL;
wjelem = NULL; wjelem = NULL;
@ -49,7 +49,8 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
bzeroflag = 1; bzeroflag = 1;
bnormflag = 0; bnormflag = 0;
quadraticflag = 0; quadraticflag = 0;
alloyflag = 0; chemflag = 0;
bnormflag = 0;
wselfallflag = 0; wselfallflag = 0;
nelements = 1; nelements = 1;
@ -101,22 +102,25 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Illegal compute snav/atom command"); error->all(FLERR,"Illegal compute snav/atom command");
quadraticflag = atoi(arg[iarg+1]); quadraticflag = atoi(arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"alloy") == 0) { } else if (strcmp(arg[iarg],"chem") == 0) {
if (iarg+2+ntypes > narg) if (iarg+2+ntypes > narg)
error->all(FLERR,"Illegal compute sna/atom command"); error->all(FLERR,"Illegal compute sna/atom command");
alloyflag = 1; chemflag = 1;
bnormflag = alloyflag;
memory->create(map,ntypes+1,"compute_sna_atom:map"); memory->create(map,ntypes+1,"compute_sna_atom:map");
nelements = force->inumeric(FLERR,arg[iarg+1]); nelements = force->inumeric(FLERR,arg[iarg+1]);
for(int i = 0; i < ntypes; i++) { for(int i = 0; i < ntypes; i++) {
int jelem = force->inumeric(FLERR,arg[iarg+2+i]); int jelem = force->inumeric(FLERR,arg[iarg+2+i]);
printf("%d %d %d %d\n",ntypes,nelements,i,jelem);
if (jelem < 0 || jelem >= nelements) if (jelem < 0 || jelem >= nelements)
error->all(FLERR,"Illegal compute snav/atom command"); error->all(FLERR,"Illegal compute snav/atom command");
map[i+1] = jelem; map[i+1] = jelem;
} }
iarg += 2+ntypes; iarg += 2+ntypes;
} else if (strcmp(arg[iarg],"wselfall") == 0) { } else if (strcmp(arg[iarg],"bnormflag") == 0) {
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command");
bnormflag = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"wselfallflag") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
error->all(FLERR,"Illegal compute snav/atom command"); error->all(FLERR,"Illegal compute snav/atom command");
wselfallflag = atoi(arg[iarg+1]); wselfallflag = atoi(arg[iarg+1]);
@ -126,7 +130,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) :
snaptr = new SNA(lmp, rfac0, twojmax, snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag, rmin0, switchflag, bzeroflag,
alloyflag, wselfallflag, nelements); chemflag, bnormflag, wselfallflag, nelements);
ncoeff = snaptr->ncoeff; ncoeff = snaptr->ncoeff;
nperdim = ncoeff; nperdim = ncoeff;
@ -236,7 +240,7 @@ void ComputeSNAVAtom::compute_peratom()
const double ztmp = x[i][2]; const double ztmp = x[i][2];
const int itype = type[i]; const int itype = type[i];
int ielem = 0; int ielem = 0;
if (alloyflag) if (chemflag)
ielem = map[itype]; ielem = map[itype];
const double radi = radelem[itype]; const double radi = radelem[itype];
@ -265,7 +269,7 @@ void ComputeSNAVAtom::compute_peratom()
const double rsq = delx*delx + dely*dely + delz*delz; const double rsq = delx*delx + dely*dely + delz*delz;
int jtype = type[j]; int jtype = type[j];
int jelem = 0; int jelem = 0;
if (alloyflag) if (chemflag)
jelem = map[jtype]; jelem = map[jtype];
if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { if (rsq < cutsq[itype][jtype]&&rsq>1e-20) {
snaptr->rij[ninside][0] = delx; snaptr->rij[ninside][0] = delx;

2
src/SNAP/compute_snav_atom.h
View File

@ -45,7 +45,7 @@ class ComputeSNAVAtom : public Compute {
double *radelem; double *radelem;
double *wjelem; double *wjelem;
int *map; // map types to [0,nelements) int *map; // map types to [0,nelements)
int nelements, alloyflag, wselfallflag; int nelements, chemflag;
class SNA* snaptr; class SNA* snaptr;
int quadraticflag; int quadraticflag;
}; };

27
src/SNAP/pair_snap.cpp
View File

@ -156,14 +156,14 @@ void PairSNAP::compute(int eflag, int vflag)
snaptr->inside[ninside] = j; snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jelem]; snaptr->wj[ninside] = wjelem[jelem];
snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
snaptr->element[ninside] = jelem; // element index for alloy snap snaptr->element[ninside] = jelem;
ninside++; ninside++;
} }
} }
// compute Ui, Yi for atom I // compute Ui, Yi for atom I
if (alloyflag) if (chemflag)
snaptr->compute_ui(ninside, ielem); snaptr->compute_ui(ninside, ielem);
else else
snaptr->compute_ui(ninside, 0); snaptr->compute_ui(ninside, 0);
@ -176,7 +176,7 @@ void PairSNAP::compute(int eflag, int vflag)
for (int jj = 0; jj < ninside; jj++) { for (int jj = 0; jj < ninside; jj++) {
int j = snaptr->inside[jj]; int j = snaptr->inside[jj];
if(alloyflag) if(chemflag)
snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj],
snaptr->rcutij[jj],jj, snaptr->element[jj]); snaptr->rcutij[jj],jj, snaptr->element[jj]);
else else
@ -328,17 +328,17 @@ void PairSNAP::compute_bispectrum()
snaptr->inside[ninside] = j; snaptr->inside[ninside] = j;
snaptr->wj[ninside] = wjelem[jelem]; snaptr->wj[ninside] = wjelem[jelem];
snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac;
snaptr->element[ninside] = jelem; // element index for alloy snap snaptr->element[ninside] = jelem;
ninside++; ninside++;
} }
} }
if (alloyflag) if (chemflag)
snaptr->compute_ui(ninside, ielem); snaptr->compute_ui(ninside, ielem);
else else
snaptr->compute_ui(ninside, 0); snaptr->compute_ui(ninside, 0);
snaptr->compute_zi(); snaptr->compute_zi();
if (alloyflag) if (chemflag)
snaptr->compute_bi(ielem); snaptr->compute_bi(ielem);
else else
snaptr->compute_bi(0); snaptr->compute_bi(0);
@ -418,7 +418,6 @@ void PairSNAP::coeff(int narg, char **arg)
ncoeffq = (ncoeff*(ncoeff+1))/2; ncoeffq = (ncoeff*(ncoeff+1))/2;
int ntmp = 1+ncoeff+ncoeffq; int ntmp = 1+ncoeff+ncoeffq;
if (ntmp != ncoeffall) { if (ntmp != ncoeffall) {
printf("ncoeffall = %d ntmp = %d ncoeff = %d \n",ncoeffall,ntmp,ncoeff);
error->all(FLERR,"Incorrect SNAP coeff file"); error->all(FLERR,"Incorrect SNAP coeff file");
} }
} }
@ -461,7 +460,7 @@ void PairSNAP::coeff(int narg, char **arg)
snaptr = new SNA(lmp, rfac0, twojmax, snaptr = new SNA(lmp, rfac0, twojmax,
rmin0, switchflag, bzeroflag, rmin0, switchflag, bzeroflag,
alloyflag, wselfallflag, nelements); chemflag, bnormflag, wselfallflag, nelements);
if (ncoeff != snaptr->ncoeff) { if (ncoeff != snaptr->ncoeff) {
if (comm->me == 0) if (comm->me == 0)
@ -653,9 +652,9 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
rmin0 = 0.0; rmin0 = 0.0;
switchflag = 1; switchflag = 1;
bzeroflag = 1; bzeroflag = 1;
bnormflag = 0;
quadraticflag = 0; quadraticflag = 0;
alloyflag = 0; chemflag = 0;
bnormflag = 0;
wselfallflag = 0; wselfallflag = 0;
chunksize = 2000; chunksize = 2000;
@ -721,8 +720,10 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
bzeroflag = atoi(keyval); bzeroflag = atoi(keyval);
else if (strcmp(keywd,"quadraticflag") == 0) else if (strcmp(keywd,"quadraticflag") == 0)
quadraticflag = atoi(keyval); quadraticflag = atoi(keyval);
else if (strcmp(keywd,"alloyflag") == 0) else if (strcmp(keywd,"chemflag") == 0)
alloyflag = atoi(keyval); chemflag = atoi(keyval);
else if (strcmp(keywd,"bnormflag") == 0)
bnormflag = atoi(keyval);
else if (strcmp(keywd,"wselfallflag") == 0) else if (strcmp(keywd,"wselfallflag") == 0)
wselfallflag = atoi(keyval); wselfallflag = atoi(keyval);
else if (strcmp(keywd,"chunksize") == 0) else if (strcmp(keywd,"chunksize") == 0)
@ -731,8 +732,6 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
error->all(FLERR,"Incorrect SNAP parameter file"); error->all(FLERR,"Incorrect SNAP parameter file");
} }
bnormflag = alloyflag;
if (rcutfacflag == 0 || twojmaxflag == 0) if (rcutfacflag == 0 || twojmaxflag == 0)
error->all(FLERR,"Incorrect SNAP parameter file"); error->all(FLERR,"Incorrect SNAP parameter file");

2
src/SNAP/pair_snap.h
View File

@ -59,7 +59,7 @@ protected:
double** bispectrum; // bispectrum components for all atoms in list double** bispectrum; // bispectrum components for all atoms in list
int *map; // mapping from atom types to elements int *map; // mapping from atom types to elements
int twojmax, switchflag, bzeroflag, bnormflag; int twojmax, switchflag, bzeroflag, bnormflag;
int alloyflag, wselfallflag; int chemflag, wselfallflag;
int chunksize; int chunksize;
double rfac0, rmin0, wj1, wj2; double rfac0, rmin0, wj1, wj2;
int rcutfacflag, twojmaxflag; // flags for required parameters int rcutfacflag, twojmaxflag; // flags for required parameters

16
src/SNAP/sna.cpp
View File

@ -109,7 +109,7 @@ using namespace MathConst;
SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in, SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in,
double rmin0_in, int switch_flag_in, int bzero_flag_in, double rmin0_in, int switch_flag_in, int bzero_flag_in,
int alloy_flag_in, int wselfall_flag_in, int nelements_in) : Pointers(lmp) int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in) : Pointers(lmp)
{ {
wself = 1.0; wself = 1.0;
@ -117,11 +117,15 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in,
rmin0 = rmin0_in; rmin0 = rmin0_in;
switch_flag = switch_flag_in; switch_flag = switch_flag_in;
bzero_flag = bzero_flag_in; bzero_flag = bzero_flag_in;
bnorm_flag = alloy_flag_in; chem_flag = chem_flag_in;
alloy_flag = alloy_flag_in; bnorm_flag = bnorm_flag_in;
wselfall_flag = wselfall_flag_in; wselfall_flag = wselfall_flag_in;
if (alloy_flag) if (bnorm_flag != chem_flag)
lmp->error->warning(FLERR, "bnormflag and chemflag are not equal."
"This is probably not what you intended");
if (chem_flag)
nelements = nelements_in; nelements = nelements_in;
else else
nelements = 1; nelements = 1;
@ -351,7 +355,7 @@ void SNA::compute_ui(int jnum, int ielem)
z0 = r / tan(theta0); z0 = r / tan(theta0);
compute_uarray(x, y, z, z0, r, j); compute_uarray(x, y, z, z0, r, j);
if (alloy_flag) if (chem_flag)
add_uarraytot(r, wj[j], rcutij[j], j, element[j]); add_uarraytot(r, wj[j], rcutij[j], j, element[j]);
else else
add_uarraytot(r, wj[j], rcutij[j], j, 0); add_uarraytot(r, wj[j], rcutij[j], j, 0);
@ -1688,7 +1692,7 @@ void SNA::compute_ncoeff()
ndoubles = nelements*nelements; ndoubles = nelements*nelements;
ntriples = nelements*nelements*nelements; ntriples = nelements*nelements*nelements;
if (alloy_flag) if (chem_flag)
ncoeff = ncount*ntriples; ncoeff = ncount*ntriples;
else else
ncoeff = ncount; ncoeff = ncount;

4
src/SNAP/sna.h
View File

@ -33,7 +33,7 @@ struct SNA_BINDICES {
class SNA : protected Pointers { class SNA : protected Pointers {
public: public:
SNA(LAMMPS*, double, int, double, int, int, int, int, int); SNA(LAMMPS*, double, int, double, int, int, int, int, int, int);
SNA(LAMMPS* lmp) : Pointers(lmp) {}; SNA(LAMMPS* lmp) : Pointers(lmp) {};
~SNA(); ~SNA();
@ -129,7 +129,7 @@ private:
int bzero_flag; // 1 if bzero subtracted from barray int bzero_flag; // 1 if bzero subtracted from barray
double* bzero; // array of B values for isolated atoms double* bzero; // array of B values for isolated atoms
int bnorm_flag; // 1 if barray divided by j+1 int bnorm_flag; // 1 if barray divided by j+1
int alloy_flag; // 1 for multi-element bispectrum components int chem_flag; // 1 for multi-element bispectrum components
int wselfall_flag; // 1 for adding wself to all element labelings int wselfall_flag; // 1 for adding wself to all element labelings
int nelements; // number of elements int nelements; // number of elements
int ndoubles; // number of multi-element pairs int ndoubles; // number of multi-element pairs