diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index 83e28b913b..9944704b77 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -30,7 +30,7 @@ Syntax * 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 * 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:: @@ -44,6 +44,15 @@ Syntax *quadraticflag* value = *0* or *1* *0* = do not 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 """""""" @@ -54,6 +63,7 @@ Examples 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 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 """"""""""" @@ -71,27 +81,26 @@ mathematical definition is given in the paper by Thompson et al. :ref:`(Thompson) ` 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) `, proposed mapping this 3D ball 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 -polar angle *theta0* defined by, +space. The radial distance *r* within *R_ii'* is mapped on to a third +polar angle :math:`\theta_0` defined by, .. 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 -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. -The natural basis for functions on the 3-sphere is formed by the 4D -hyperspherical harmonics *U\^j_m,m'(theta, phi, theta0).* These -functions are better known as *D\^j_m,m',* the elements of the Wigner +The natural basis for functions on the 3-sphere is formed by the +representatives of *SU(2)*, the matrices :math:`U^j_{m,m'}(\theta, \phi, \theta_0)`. +These functions are better known as :math:`D^j_{m,m'}`, the elements of the Wigner *D*\ -matrices :ref:`(Meremianin `, -:ref:`Varshalovich) `. - +:ref:`Varshalovich `, :ref:`Mason) ` 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 radial distance. Expanding this density function as a generalized @@ -100,20 +109,20 @@ coefficient as .. 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 -to distinguish atoms of different types, while the central atom is -arbitrarily assigned a unit weight. The function *fc(r)* ensures that +The :math:`w_{\mu_{i'}}` neighbor weights are dimensionless numbers that depend on +:math:`\mu_{i'}`, the SNAP element of atom *i'*, while the central atom is +arbitrarily assigned a unit weight. The function :math:`f_c(r)` ensures that the contribution of each neighbor atom goes smoothly to zero at -*R_ii'*: +:math:`R_{ii'}`: .. math:: 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'} -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 under rotation of the polar coordinate frame. However, the following scalar triple products of expansion coefficients can be shown to be @@ -128,7 +137,8 @@ real-valued and invariant under rotation :ref:`(Bartok) `. {j_2} {m_2} {m'_2} \end{array}} 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 2-sphere. These invariants are the components of the bispectrum and 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 characterizes the correlations of two points. The lowest-order components describe the coarsest features of the density function, -while higher-order components reflect finer detail. Note that the -central atom is included in the expansion, so three point-correlations -can be either due to three neighbors, or two neighbors and the central -atom. +while higher-order components reflect finer detail. Each bispectrum +component contains terms that depend on the positions of up to 4 +atoms (3 neighbors and the central atom). Compute *snad/atom* calculates the derivative of the bispectrum components -summed separately for each atom type: +summed separately for each LAMMPS atom type: .. 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 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*\ , 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 @@ -201,8 +210,8 @@ The argument *rcutfac* is a scale factor that controls the ratio of atomic radius to radial cutoff distance. The argument *rfac0* and the optional keyword *rmin0* define the -linear mapping from radial distance to polar angle *theta0* on the -3-sphere. +linear mapping from radial distance to polar angle :math:`theta_0` on the +3-sphere, given above. The argument *twojmax* defines which 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. 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 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*\ . 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 of bispectrum components with itself. See section below on output for a detailed explanation of the number of quadratic terms and the 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. ` +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:: 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): 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:: 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 arranged into *ntypes* blocks, listed in order of atom type *I*\ . Each block contains three sub-blocks corresponding to the *x*\ , *y*\ , and *z* -components of the atom position. Each of these sub-blocks contains -one column for each bispectrum component, the same as for compute -*sna/atom* +components of the atom position. Each of these sub-blocks contains *K* +columns for the *K* bispectrum components, the same as for compute *sna/atom* Compute *snav/atom* evaluates a per-atom array. The columns are arranged into *ntypes* blocks, listed in order of atom type *I*\ . Each block contains six sub-blocks corresponding to the *xx*\ , *yy*\ , *zz*\ , *yz*\ , *xz*\ , and *xy* components of the virial tensor in Voigt -notation. Each of these sub-blocks contains one column for each -bispectrum component, the same as for compute *sna/atom* +notation. Each of these sub-blocks contains *K* +columns for the *K* bispectrum components, the same as for compute *sna/atom* Compute *snap* evaluates a global array. 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, 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 from a compute as input. See the :doc:`Howto output ` doc 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 -LAMMPS was built with that package. See the :doc:`Build package ` doc page for more info. +LAMMPS was built with that package. See the :doc:`Build package ` +doc page for more info. Related commands """""""""""""""" @@ -332,6 +410,7 @@ Default The optional keyword defaults are *rmin0* = 0, *switchflag* = 1, *bzeroflag* = 1, *quadraticflag* = 0, +*bnormflag* = 0, *wselfallflag* = 0 ---------- @@ -352,3 +431,12 @@ available at `arXiv:1409.3880 `_ **(Varshalovich)** Varshalovich, Moskalev, Khersonskii, Quantum Theory 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) diff --git a/doc/src/pair_snap.rst b/doc/src/pair_snap.rst index 0770b5a675..a46df30609 100644 --- a/doc/src/pair_snap.rst +++ b/doc/src/pair_snap.rst @@ -24,27 +24,30 @@ Examples Description """"""""""" -Pair style *snap* computes interactions using the spectral -neighbor analysis potential (SNAP) :ref:`(Thompson) `. -Like the GAP framework of Bartok et al. :ref:`(Bartok2010) `, -:ref:`(Bartok2013) ` which uses bispectrum components +Pair style *snap* defines the spectral +neighbor analysis potential (SNAP), a machine-learning +interatomic potential :ref:`(Thompson) `. +Like the GAP framework of Bartok et al. :ref:`(Bartok2010) `, +SNAP uses bispectrum components to characterize the local neighborhood of each atom in a very general way. The mathematical definition of the -bispectrum calculation used by SNAP is identical -to that used by :doc:`compute sna/atom `. +bispectrum calculation and its derivatives w.r.t. atom positions +is identical to that used by :doc:`compute snap `, +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 atom energies. The energy of atom *i* is expressed as a weighted sum over bispectrum components. .. 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*\ , -and :math:`\beta_k^{\alpha_i}` is the corresponding linear coefficient -that depends on :math:\alpha_i`, the SNAP element of atom *i*\ . The +and :math:`\beta_k^{\mu_i}` is the corresponding linear coefficient +that depends on :math:`\mu_i`, the SNAP element of atom *i*\ . The 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. The bispectrum calculation is described in more detail in :doc:`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 keyword/value pair. The required keywords are *rcutfac* and *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 * *rfac0* = 0.99363 * *rmin0* = 0.0 -* *switchflag* = 0 +* *switchflag* = 1 * *bzeroflag* = 1 -* *quadraticflag* = 1 +* *quadraticflag* = 0 +* *chemflag* = 0 +* *bnormflag* = 0 +* *wselfallflag* = 0 * *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) `. + +.. 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) `. + +.. 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 pair style *snap* with the KOKKOS package and is ignored otherwise. This keyword controls @@ -159,10 +196,6 @@ into two passes. Detailed definitions for all the other keywords are given on the :doc:`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:: The previously used *diagonalstyle* keyword was removed in 2019, @@ -221,7 +254,8 @@ Related commands :doc:`compute sna/atom `, :doc:`compute snad/atom `, -:doc:`compute snav/atom ` +:doc:`compute snav/atom `, +:doc:`compute snap ` **Default:** none @@ -235,6 +269,10 @@ Related commands **(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) diff --git a/examples/snap/InP_Cusentino_PRB2020.snap b/examples/snap/InP_Cusentino_PRB2020.snap deleted file mode 100644 index bada0b04a9..0000000000 --- a/examples/snap/InP_Cusentino_PRB2020.snap +++ /dev/null @@ -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 diff --git a/examples/snap/InP_Cusentino_PRB2020.snapparam b/examples/snap/InP_Cusentino_PRB2020.snapparam deleted file mode 100644 index 521b8a4521..0000000000 --- a/examples/snap/InP_Cusentino_PRB2020.snapparam +++ /dev/null @@ -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 - \ No newline at end of file diff --git a/examples/snap/InP_JCPA2020.snap b/examples/snap/InP_JCPA2020.snap new file mode 120000 index 0000000000..d4b8ea6baa --- /dev/null +++ b/examples/snap/InP_JCPA2020.snap @@ -0,0 +1 @@ +../../potentials/InP_JCPA2020.snap \ No newline at end of file diff --git a/examples/snap/InP_JCPA2020.snapcoeff b/examples/snap/InP_JCPA2020.snapcoeff new file mode 120000 index 0000000000..c99cb558b1 --- /dev/null +++ b/examples/snap/InP_JCPA2020.snapcoeff @@ -0,0 +1 @@ +../../potentials/InP_JCPA2020.snapcoeff \ No newline at end of file diff --git a/examples/snap/InP_JCPA2020.snapparam b/examples/snap/InP_JCPA2020.snapparam new file mode 120000 index 0000000000..1c090c8001 --- /dev/null +++ b/examples/snap/InP_JCPA2020.snapparam @@ -0,0 +1 @@ +../../potentials/InP_JCPA2020.snapparam \ No newline at end of file diff --git a/examples/snap/in.snap.InP.PRB2020 b/examples/snap/in.snap.InP.JCPA2020 similarity index 92% rename from examples/snap/in.snap.InP.PRB2020 rename to examples/snap/in.snap.InP.JCPA2020 index a4a9bfc7bf..dd16301a01 100644 --- a/examples/snap/in.snap.InP.PRB2020 +++ b/examples/snap/in.snap.InP.JCPA2020 @@ -6,7 +6,6 @@ variable nsteps index 100 variable nrep equal 4 variable a equal 5.83 units metal -atom_style charge # generate the box and atom positions using a FCC lattice @@ -26,7 +25,7 @@ mass 2 30.98 # choose potential -include InP_Cusentino_PRB2020.snap +include InP_JCPA2020.snap # Setup output diff --git a/examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.1 b/examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.1 new file mode 100644 index 0000000000..094e039420 --- /dev/null +++ b/examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.1 @@ -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 diff --git a/examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.4 b/examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.4 new file mode 100644 index 0000000000..109ed315d2 --- /dev/null +++ b/examples/snap/log.01Jun20.snap.InP_JCPA2020.g++.4 @@ -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 diff --git a/potentials/InP_JCPA2020.snap b/potentials/InP_JCPA2020.snap new file mode 100644 index 0000000000..6c443df739 --- /dev/null +++ b/potentials/InP_JCPA2020.snap @@ -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 diff --git a/examples/snap/InP_Cusentino_PRB2020.snapcoeff b/potentials/InP_JCPA2020.snapcoeff similarity index 98% rename from examples/snap/InP_Cusentino_PRB2020.snapcoeff rename to potentials/InP_JCPA2020.snapcoeff index a589546cf9..35f0ab8d42 100644 --- a/examples/snap/InP_Cusentino_PRB2020.snapcoeff +++ b/potentials/InP_JCPA2020.snapcoeff @@ -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 In 3.81205 1 diff --git a/potentials/InP_JCPA2020.snapparam b/potentials/InP_JCPA2020.snapparam new file mode 100644 index 0000000000..83c576f8f4 --- /dev/null +++ b/potentials/InP_JCPA2020.snapparam @@ -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 + \ No newline at end of file diff --git a/src/KOKKOS/pair_snap_kokkos_impl.h b/src/KOKKOS/pair_snap_kokkos_impl.h index e992396f7b..5c2eb83746 100644 --- a/src/KOKKOS/pair_snap_kokkos_impl.h +++ b/src/KOKKOS/pair_snap_kokkos_impl.h @@ -515,7 +515,7 @@ void PairSNAPKokkos::coeff(int narg, char **arg) Kokkos::deep_copy(d_map,h_map); snaKK = SNAKokkos(rfac0,twojmax, - rmin0,switchflag,bzeroflag,alloyflag,wselfallflag,nelements); + rmin0,switchflag,bzeroflag,chemflag,bnormflag,wselfallflag,nelements); snaKK.grow_rij(0,0); snaKK.init(); } @@ -584,7 +584,7 @@ void PairSNAPKokkos::operator() (TagPairSNAPComputeNeigh,const typen my_sna.inside(ii,offset) = j; my_sna.wj(ii,offset) = d_wjelem[elem_j]; my_sna.rcutij(ii,offset) = (radi + d_radelem[elem_j])*rcutfac; - if (alloyflag) + if (chemflag) my_sna.element(ii,offset) = elem_j; else my_sna.element(ii,offset) = 0; diff --git a/src/KOKKOS/sna_kokkos.h b/src/KOKKOS/sna_kokkos.h index f72fa6df37..f9dbb07c92 100644 --- a/src/KOKKOS/sna_kokkos.h +++ b/src/KOKKOS/sna_kokkos.h @@ -210,7 +210,7 @@ inline int switch_flag; // Chem snap flags - int alloy_flag; + int chem_flag; int bnorm_flag; int nelements; int ndoubles; diff --git a/src/KOKKOS/sna_kokkos_impl.h b/src/KOKKOS/sna_kokkos_impl.h index 72185c7fb9..0250c43058 100644 --- a/src/KOKKOS/sna_kokkos_impl.h +++ b/src/KOKKOS/sna_kokkos_impl.h @@ -29,7 +29,7 @@ template inline SNAKokkos::SNAKokkos(double rfac0_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; @@ -38,13 +38,13 @@ SNAKokkos::SNAKokkos(double rfac0_in, switch_flag = switch_flag_in; bzero_flag = bzero_flag_in; - alloy_flag = alloy_flag_in; - wselfall_flag = wselfall_flag_in; - if (alloy_flag) + chem_flag = chem_flag_in; + if (chem_flag) nelements = nelements_in; else nelements = 1; - bnorm_flag = alloy_flag_in; + bnorm_flag = bnorm_flag_in; + wselfall_flag = wselfall_flag_in; twojmax = twojmax_in; @@ -283,7 +283,7 @@ void SNAKokkos::pre_ui(const typename Kokkos::TeamPolicy // if m is on the "diagonal", initialize it with the self energy. // Otherwise zero it out 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; }); @@ -1620,7 +1620,7 @@ int SNAKokkos::compute_ncoeff() ndoubles = nelements*nelements; ntriples = nelements*nelements*nelements; - if (alloy_flag) ncount *= ntriples; + if (chem_flag) ncount *= ntriples; return ncount; } diff --git a/src/SNAP/compute_sna_atom.cpp b/src/SNAP/compute_sna_atom.cpp index 0d81776037..a90bfbfa28 100644 --- a/src/SNAP/compute_sna_atom.cpp +++ b/src/SNAP/compute_sna_atom.cpp @@ -34,7 +34,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : radelem(NULL), wjelem(NULL) { double rmin0, rfac0; - int twojmax, switchflag, bzeroflag, bnormflag; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; radelem = NULL; wjelem = NULL; @@ -50,7 +50,8 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : bzeroflag = 1; bnormflag = 0; quadraticflag = 0; - alloyflag = 0; + chemflag = 0; + bnormflag = 0; wselfallflag = 0; nelements = 1; @@ -108,22 +109,25 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute sna/atom command"); quadraticflag = atoi(arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"alloy") == 0) { + } else if (strcmp(arg[iarg],"chem") == 0) { if (iarg+2+ntypes > narg) error->all(FLERR,"Illegal compute sna/atom command"); - alloyflag = 1; - bnormflag = alloyflag; + chemflag = 1; memory->create(map,ntypes+1,"compute_sna_atom:map"); nelements = force->inumeric(FLERR,arg[iarg+1]); for(int i = 0; i < ntypes; 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) error->all(FLERR,"Illegal compute sna/atom command"); map[i+1] = jelem; } 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) error->all(FLERR,"Illegal compute sna/atom command"); wselfallflag = atoi(arg[iarg+1]); @@ -133,7 +137,7 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - alloyflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, nelements); ncoeff = snaptr->ncoeff; size_peratom_cols = ncoeff; @@ -229,7 +233,7 @@ void ComputeSNAAtom::compute_peratom() const double ztmp = x[i][2]; const int itype = type[i]; int ielem = 0; - if (alloyflag) + if (chemflag) ielem = map[itype]; const double radi = radelem[itype]; const int* const jlist = firstneigh[i]; @@ -254,7 +258,7 @@ void ComputeSNAAtom::compute_peratom() const double rsq = delx*delx + dely*dely + delz*delz; int jtype = type[j]; int jelem = 0; - if (alloyflag) + if (chemflag) int jelem = map[jtype]; if (rsq < cutsq[itype][jtype] && rsq>1e-20) { snaptr->rij[ninside][0] = delx; diff --git a/src/SNAP/compute_sna_atom.h b/src/SNAP/compute_sna_atom.h index 870e7c0744..1ad325729d 100644 --- a/src/SNAP/compute_sna_atom.h +++ b/src/SNAP/compute_sna_atom.h @@ -43,7 +43,7 @@ class ComputeSNAAtom : public Compute { double *radelem; double *wjelem; int * map; // map types to [0,nelements) - int nelements, alloyflag, wselfallflag; + int nelements, chemflag; class SNA* snaptr; double cutmax; int quadraticflag; diff --git a/src/SNAP/compute_snad_atom.cpp b/src/SNAP/compute_snad_atom.cpp index d500955384..cd0c0cbc41 100644 --- a/src/SNAP/compute_snad_atom.cpp +++ b/src/SNAP/compute_snad_atom.cpp @@ -34,7 +34,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : radelem(NULL), wjelem(NULL) { double rfac0, rmin0; - int twojmax, switchflag, bzeroflag, bnormflag; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; radelem = NULL; wjelem = NULL; @@ -50,7 +50,8 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : bzeroflag = 1; bnormflag = 0; quadraticflag = 0; - alloyflag = 0; + chemflag = 0; + bnormflag = 0; wselfallflag = 0; nelements = 1; @@ -106,22 +107,25 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snad/atom command"); quadraticflag = atoi(arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"alloy") == 0) { + } else if (strcmp(arg[iarg],"chem") == 0) { if (iarg+2+ntypes > narg) error->all(FLERR,"Illegal compute snad/atom command"); - alloyflag = 1; - bnormflag = alloyflag; + chemflag = 1; memory->create(map,ntypes+1,"compute_snad_atom:map"); nelements = force->inumeric(FLERR,arg[iarg+1]); for(int i = 0; i < ntypes; 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) error->all(FLERR,"Illegal compute snad/atom command"); map[i+1] = jelem; } 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) error->all(FLERR,"Illegal compute snad/atom command"); wselfallflag = atoi(arg[iarg+1]); @@ -131,7 +135,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - alloyflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, nelements); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -241,7 +245,7 @@ void ComputeSNADAtom::compute_peratom() const double ztmp = x[i][2]; const int itype = type[i]; int ielem = 0; - if (alloyflag) + if (chemflag) ielem = map[itype]; const double radi = radelem[itype]; const int* const jlist = firstneigh[i]; @@ -272,7 +276,7 @@ void ComputeSNADAtom::compute_peratom() const double rsq = delx*delx + dely*dely + delz*delz; int jtype = type[j]; int jelem = 0; - if (alloyflag) + if (chemflag) jelem = map[jtype]; if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { snaptr->rij[ninside][0] = delx; diff --git a/src/SNAP/compute_snad_atom.h b/src/SNAP/compute_snad_atom.h index 6a1155adac..15609bdfdc 100644 --- a/src/SNAP/compute_snad_atom.h +++ b/src/SNAP/compute_snad_atom.h @@ -45,7 +45,7 @@ class ComputeSNADAtom : public Compute { double *radelem; double *wjelem; int *map; // map types to [0,nelements) - int nelements, alloyflag, wselfallflag; + int nelements, chemflag; class SNA* snaptr; double cutmax; int quadraticflag; diff --git a/src/SNAP/compute_snap.cpp b/src/SNAP/compute_snap.cpp index 251a823799..c9e1006571 100644 --- a/src/SNAP/compute_snap.cpp +++ b/src/SNAP/compute_snap.cpp @@ -41,7 +41,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : extarray = 0; double rfac0, rmin0; - int twojmax, switchflag, bzeroflag; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; radelem = NULL; wjelem = NULL; @@ -56,7 +56,8 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : switchflag = 1; bzeroflag = 1; quadraticflag = 0; - alloyflag = 0; + chemflag = 0; + bnormflag = 0; wselfallflag = 0; nelements = 1; @@ -112,20 +113,24 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snap command"); quadraticflag = atoi(arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"alloyflag") == 0) { + } else if (strcmp(arg[iarg],"chem") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute snap command"); - alloyflag = 1; + chemflag = 1; memory->create(map,ntypes+1,"compute_snap:map"); nelements = force->inumeric(FLERR,arg[iarg+1]); for(int i = 0; i < ntypes; 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) error->all(FLERR,"Illegal compute snap command"); map[i+1] = jelem; } 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) { if (iarg+2 > narg) 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, rmin0, switchflag, bzeroflag, - alloyflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, nelements); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -168,7 +173,7 @@ ComputeSnap::~ComputeSnap() memory->destroy(cutsq); 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 int itype = type[i]; int ielem = 0; - if (alloyflag) + if (chemflag) ielem = map[itype]; const double radi = radelem[itype]; const int* const jlist = firstneigh[i]; @@ -328,7 +333,7 @@ void ComputeSnap::compute_array() const double rsq = delx*delx + dely*dely + delz*delz; int jtype = type[j]; int jelem = 0; - if (alloyflag) + if (chemflag) jelem = map[jtype]; if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { snaptr->rij[ninside][0] = delx; diff --git a/src/SNAP/compute_snap.h b/src/SNAP/compute_snap.h index 4e8219c693..7151839c8c 100644 --- a/src/SNAP/compute_snap.h +++ b/src/SNAP/compute_snap.h @@ -45,7 +45,7 @@ class ComputeSnap : public Compute { double *radelem; double *wjelem; int *map; // map types to [0,nelements) - int nelements, alloyflag, wselfallflag; + int nelements, chemflag; class SNA* snaptr; double cutmax; int quadraticflag; diff --git a/src/SNAP/compute_snav_atom.cpp b/src/SNAP/compute_snav_atom.cpp index ffec79fe95..2d33c35563 100644 --- a/src/SNAP/compute_snav_atom.cpp +++ b/src/SNAP/compute_snav_atom.cpp @@ -33,7 +33,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : radelem(NULL), wjelem(NULL) { double rfac0, rmin0; - int twojmax, switchflag, bzeroflag, bnormflag; + int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; radelem = NULL; wjelem = NULL; @@ -49,7 +49,8 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : bzeroflag = 1; bnormflag = 0; quadraticflag = 0; - alloyflag = 0; + chemflag = 0; + bnormflag = 0; wselfallflag = 0; nelements = 1; @@ -101,22 +102,25 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Illegal compute snav/atom command"); quadraticflag = atoi(arg[iarg+1]); iarg += 2; - } else if (strcmp(arg[iarg],"alloy") == 0) { + } else if (strcmp(arg[iarg],"chem") == 0) { if (iarg+2+ntypes > narg) error->all(FLERR,"Illegal compute sna/atom command"); - alloyflag = 1; - bnormflag = alloyflag; + chemflag = 1; memory->create(map,ntypes+1,"compute_sna_atom:map"); nelements = force->inumeric(FLERR,arg[iarg+1]); for(int i = 0; i < ntypes; 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) error->all(FLERR,"Illegal compute snav/atom command"); map[i+1] = jelem; } 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) error->all(FLERR,"Illegal compute snav/atom command"); wselfallflag = atoi(arg[iarg+1]); @@ -126,7 +130,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - alloyflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, nelements); ncoeff = snaptr->ncoeff; nperdim = ncoeff; @@ -236,7 +240,7 @@ void ComputeSNAVAtom::compute_peratom() const double ztmp = x[i][2]; const int itype = type[i]; int ielem = 0; - if (alloyflag) + if (chemflag) ielem = map[itype]; const double radi = radelem[itype]; @@ -265,7 +269,7 @@ void ComputeSNAVAtom::compute_peratom() const double rsq = delx*delx + dely*dely + delz*delz; int jtype = type[j]; int jelem = 0; - if (alloyflag) + if (chemflag) jelem = map[jtype]; if (rsq < cutsq[itype][jtype]&&rsq>1e-20) { snaptr->rij[ninside][0] = delx; diff --git a/src/SNAP/compute_snav_atom.h b/src/SNAP/compute_snav_atom.h index f0789282d5..857671388b 100644 --- a/src/SNAP/compute_snav_atom.h +++ b/src/SNAP/compute_snav_atom.h @@ -45,7 +45,7 @@ class ComputeSNAVAtom : public Compute { double *radelem; double *wjelem; int *map; // map types to [0,nelements) - int nelements, alloyflag, wselfallflag; + int nelements, chemflag; class SNA* snaptr; int quadraticflag; }; diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp index f25743176e..63867f1691 100644 --- a/src/SNAP/pair_snap.cpp +++ b/src/SNAP/pair_snap.cpp @@ -156,14 +156,14 @@ void PairSNAP::compute(int eflag, int vflag) snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; - snaptr->element[ninside] = jelem; // element index for alloy snap + snaptr->element[ninside] = jelem; ninside++; } } // compute Ui, Yi for atom I - if (alloyflag) + if (chemflag) snaptr->compute_ui(ninside, ielem); else snaptr->compute_ui(ninside, 0); @@ -176,7 +176,7 @@ void PairSNAP::compute(int eflag, int vflag) for (int jj = 0; jj < ninside; jj++) { int j = snaptr->inside[jj]; - if(alloyflag) + if(chemflag) snaptr->compute_duidrj(snaptr->rij[jj], snaptr->wj[jj], snaptr->rcutij[jj],jj, snaptr->element[jj]); else @@ -328,17 +328,17 @@ void PairSNAP::compute_bispectrum() snaptr->inside[ninside] = j; snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; - snaptr->element[ninside] = jelem; // element index for alloy snap + snaptr->element[ninside] = jelem; ninside++; } } - if (alloyflag) + if (chemflag) snaptr->compute_ui(ninside, ielem); else snaptr->compute_ui(ninside, 0); snaptr->compute_zi(); - if (alloyflag) + if (chemflag) snaptr->compute_bi(ielem); else snaptr->compute_bi(0); @@ -418,7 +418,6 @@ void PairSNAP::coeff(int narg, char **arg) ncoeffq = (ncoeff*(ncoeff+1))/2; int ntmp = 1+ncoeff+ncoeffq; if (ntmp != ncoeffall) { - printf("ncoeffall = %d ntmp = %d ncoeff = %d \n",ncoeffall,ntmp,ncoeff); error->all(FLERR,"Incorrect SNAP coeff file"); } } @@ -461,7 +460,7 @@ void PairSNAP::coeff(int narg, char **arg) snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, - alloyflag, wselfallflag, nelements); + chemflag, bnormflag, wselfallflag, nelements); if (ncoeff != snaptr->ncoeff) { if (comm->me == 0) @@ -653,9 +652,9 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) rmin0 = 0.0; switchflag = 1; bzeroflag = 1; - bnormflag = 0; quadraticflag = 0; - alloyflag = 0; + chemflag = 0; + bnormflag = 0; wselfallflag = 0; chunksize = 2000; @@ -721,8 +720,10 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) bzeroflag = atoi(keyval); else if (strcmp(keywd,"quadraticflag") == 0) quadraticflag = atoi(keyval); - else if (strcmp(keywd,"alloyflag") == 0) - alloyflag = atoi(keyval); + else if (strcmp(keywd,"chemflag") == 0) + chemflag = atoi(keyval); + else if (strcmp(keywd,"bnormflag") == 0) + bnormflag = atoi(keyval); else if (strcmp(keywd,"wselfallflag") == 0) wselfallflag = atoi(keyval); 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"); } - bnormflag = alloyflag; - if (rcutfacflag == 0 || twojmaxflag == 0) error->all(FLERR,"Incorrect SNAP parameter file"); diff --git a/src/SNAP/pair_snap.h b/src/SNAP/pair_snap.h index 3ed1c1ff9d..d60000da3a 100644 --- a/src/SNAP/pair_snap.h +++ b/src/SNAP/pair_snap.h @@ -59,7 +59,7 @@ protected: double** bispectrum; // bispectrum components for all atoms in list int *map; // mapping from atom types to elements int twojmax, switchflag, bzeroflag, bnormflag; - int alloyflag, wselfallflag; + int chemflag, wselfallflag; int chunksize; double rfac0, rmin0, wj1, wj2; int rcutfacflag, twojmaxflag; // flags for required parameters diff --git a/src/SNAP/sna.cpp b/src/SNAP/sna.cpp index cbd05645a5..e003ec7ca8 100644 --- a/src/SNAP/sna.cpp +++ b/src/SNAP/sna.cpp @@ -109,7 +109,7 @@ using namespace MathConst; SNA::SNA(LAMMPS* lmp, double rfac0_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) : Pointers(lmp) + int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in) : Pointers(lmp) { wself = 1.0; @@ -117,11 +117,15 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in, rmin0 = rmin0_in; switch_flag = switch_flag_in; bzero_flag = bzero_flag_in; - bnorm_flag = alloy_flag_in; - alloy_flag = alloy_flag_in; + chem_flag = chem_flag_in; + bnorm_flag = bnorm_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; else nelements = 1; @@ -351,7 +355,7 @@ void SNA::compute_ui(int jnum, int ielem) z0 = r / tan(theta0); compute_uarray(x, y, z, z0, r, j); - if (alloy_flag) + if (chem_flag) add_uarraytot(r, wj[j], rcutij[j], j, element[j]); else add_uarraytot(r, wj[j], rcutij[j], j, 0); @@ -1688,7 +1692,7 @@ void SNA::compute_ncoeff() ndoubles = nelements*nelements; ntriples = nelements*nelements*nelements; - if (alloy_flag) + if (chem_flag) ncoeff = ncount*ntriples; else ncoeff = ncount; diff --git a/src/SNAP/sna.h b/src/SNAP/sna.h index 873860c80b..30c4aa6bbc 100644 --- a/src/SNAP/sna.h +++ b/src/SNAP/sna.h @@ -33,7 +33,7 @@ struct SNA_BINDICES { class SNA : protected Pointers { 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(); @@ -129,7 +129,7 @@ private: int bzero_flag; // 1 if bzero subtracted from barray double* bzero; // array of B values for isolated atoms 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 nelements; // number of elements int ndoubles; // number of multi-element pairs