diff --git a/doc/src/compute_sna_atom.rst b/doc/src/compute_sna_atom.rst index e43e61cd00..54a6df02a2 100644 --- a/doc/src/compute_sna_atom.rst +++ b/doc/src/compute_sna_atom.rst @@ -33,7 +33,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* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* +* keyword = *rmin0* or *switchflag* or *bzeroflag* or *quadraticflag* or *chem* or *bnormflag* or *wselfallflag* or *bikflag* or *switchinnerflag* or *sinner* or *dinner* .. parsed-literal:: @@ -59,9 +59,13 @@ Syntax *bikflag* value = *0* or *1* (only implemented for compute snap) *0* = per-atom bispectrum descriptors are summed over atoms *1* = per-atom bispectrum descriptors are not summed over atoms - *switchinnerflag* values = *rinnerlist* *drinnerlist* - *rinnerlist* = *ntypes* values of rinner (distance units) - *drinnerlist* = *ntypes* values of drinner (distance units) + *switchinnerflag* value = *0* or *1* + *0* = do not use inner switching function + *1* = use inner switching function + *sinner* values = *sinnerlist* + *sinnerlist* = *ntypes* values of *Sinner* (distance units) + *dinner* values = *dinnerlist* + *dinnerlist* = *ntypes* values of *Dinner* (distance units) Examples """""""" @@ -73,7 +77,7 @@ Examples 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 - compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1.1 1.3 0.5 0.6 + compute snap all snap 1.0 0.99363 6 3.81 3.83 1.0 0.93 switchinnerflag 1 sinner 1.35 1.6 dinner 0.25 0.3 Description """"""""""" @@ -312,25 +316,29 @@ the resulting bispectrum rows are :math:`B_{i,k}` instead of just :math:`B_k`. In this case, the entries in the final column for these rows are set to zero. -The keyword *switchinnerflag* activates an additional radial switching +The keyword *switchinnerflag* with value 1 +activates an additional radial switching function similar to :math:`f_c(r)` above, but acting to switch off smoothly contributions from neighbor atoms at short separation distances. This is useful when SNAP is used in combination with a simple -repulsive potential. The keyword is followed by the *ntypes* -values for :math:`r_{inner}` and the *ntypes* -values for :math:`\Delta r_{inner}`. For a neighbor atom at +repulsive potential. For a neighbor atom at distance :math:`r`, its contribution is scaled by a multiplicative factor :math:`f_{inner}(r)` defined as follows: .. math:: - = & 0, r \leq r_{inner} \\ - f_{inner}(r) = & \frac{1}{2}(1 - \cos(\pi \frac{r-r_{inner}}{\Delta r_{inner}})), r_{inner} < r \leq r_{inner} + \Delta r_{inner} \\ - = & 1, r > r_{inner} + \Delta r_{inner} + = & 0, r \leq S_{inner} - D_{inner} \\ + f_{inner}(r) = & \frac{1}{2}(1 - \cos(\frac{\pi}{2} (1 + \frac{r-S_{inner}}{D_{inner}})), S_{inner} - D_{inner} < r \leq S_{inner} + D_{inner} \\ + = & 1, r > S_{inner} + D_{inner} -The values of :math:`r_{inner}` and :math:`\Delta r_{inner}` are -the arithmetic means of the values for the central atom of type I -and the neighbor atom of type J. +where the switching region is centered at :math:`S_{inner}` and it extends a distance :math:`D_{inner}` +to the left and to the right of this. +With this option, additional keywords *sinner* and *dinner* must be used, +each followed by *ntypes* +values for :math:`S_{inner}` and :math:`D_{inner}`, respectively. +When the central atom and the neighbor atom have different types, +the values of :math:`S_{inner}` and :math:`D_{inner}` are +the arithmetic means of the values for both types. .. note:: @@ -450,7 +458,7 @@ Default The optional keyword defaults are *rmin0* = 0, *switchflag* = 1, *bzeroflag* = 1, *quadraticflag* = 0, -*bnormflag* = 0, *wselfallflag* = 0 +*bnormflag* = 0, *wselfallflag* = 0, *switchinnerflag* = 0, ---------- diff --git a/doc/src/pair_snap.rst b/doc/src/pair_snap.rst index 350561fe4d..ebedb288c1 100644 --- a/doc/src/pair_snap.rst +++ b/doc/src/pair_snap.rst @@ -136,7 +136,7 @@ keyword/value pair. The required keywords are *rcutfac* and *twojmax*\ . Optional keywords are *rfac0*, *rmin0*, *switchflag*, *bzeroflag*, *quadraticflag*, *chemflag*, *bnormflag*, *wselfallflag*, *switchinnerflag*, -*rinner*, *drinner*, *chunksize*, and *parallelthresh*\ . +*sinner*, *dinner*, *chunksize*, and *parallelthresh*\ . The default values for these keywords are @@ -152,6 +152,9 @@ The default values for these keywords are * *chunksize* = 32768 * *parallelthresh* = 8192 +For detailed definitions of all of these keywords, +see the :doc:`compute sna/atom ` doc page. + 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 @@ -194,7 +197,7 @@ pair_coeff command, to avoid ambiguity in the number of coefficients. The keyword *switchinnerflag* activates an additional switching function that smoothly turns off contributions to the SNAP potential from neighbor atoms at short separations. If *switchinnerflag* is set to 1 then -the additional keywords *rinner* and *drinner* must also be provided. +the additional keywords *sinner* and *dinner* must also be provided. Each of these is followed by *nelements* values, where *nelements* is the number of unique elements appearing in appearing in the LAMMPS pair_coeff command. The element order should correspond to the order @@ -217,9 +220,6 @@ already large enough to saturate the GPU threads. Extra parallelism will be performed if the *chunksize* (or total number of atoms per GPU) is smaller than *parallelthresh*. -Detailed definitions for all the other keywords -are given on the :doc:`compute sna/atom ` doc page. - .. note:: The previously used *diagonalstyle* keyword was removed in 2019, diff --git a/src/KOKKOS/pair_snap_kokkos.h b/src/KOKKOS/pair_snap_kokkos.h index bdad113c18..2cd0003478 100644 --- a/src/KOKKOS/pair_snap_kokkos.h +++ b/src/KOKKOS/pair_snap_kokkos.h @@ -73,7 +73,7 @@ struct TagPairSNAPComputeDeidrjCPU{}; template class PairSNAPKokkos : public PairSNAP { -public: + public: enum {EnabledNeighFlags=FULL|HALF|HALFTHREAD}; enum {COUL_FLAG=0}; typedef DeviceType device_type; @@ -228,7 +228,7 @@ public: const F_FLOAT &fx, const F_FLOAT &fy, const F_FLOAT &fz, const F_FLOAT &delx, const F_FLOAT &dely, const F_FLOAT &delz) const; -protected: + protected: typename AT::t_neighbors_2d d_neighbors; typename AT::t_int_1d_randomread d_ilist; typename AT::t_int_1d_randomread d_numneigh; @@ -252,9 +252,9 @@ protected: void allocate() override; //void read_files(char *, char *); /*template -inline int equal(double* x,double* y); + inline int equal(double* x,double* y); template -inline double dist2(double* x,double* y); + inline double dist2(double* x,double* y); double extra_cutoff(); void load_balance(); void set_sna_to_shared(int snaid,int i); @@ -280,6 +280,8 @@ inline double dist2(double* x,double* y); Kokkos::View d_radelem; // element radii Kokkos::View d_wjelem; // elements weights Kokkos::View d_coeffelem; // element bispectrum coefficients + Kokkos::View d_sinnerelem; // element inner cutoff midpoint + Kokkos::View d_dinnerelem; // element inner cutoff half-width Kokkos::View d_map; // mapping from atom types to elements Kokkos::View d_ninside; // ninside for all atoms in list Kokkos::View d_beta; // betas for all atoms in list @@ -329,10 +331,10 @@ inline double dist2(double* x,double* y); template class PairSNAPKokkosDevice : public PairSNAPKokkos { -private: + private: using Base = PairSNAPKokkos; -public: + public: PairSNAPKokkosDevice(class LAMMPS *); @@ -348,10 +350,10 @@ public: template class PairSNAPKokkosHost : public PairSNAPKokkos { -private: + private: using Base = PairSNAPKokkos; -public: + public: PairSNAPKokkosHost(class LAMMPS *); diff --git a/src/KOKKOS/pair_snap_kokkos_impl.h b/src/KOKKOS/pair_snap_kokkos_impl.h index 6cf40b31a5..74264ca47e 100644 --- a/src/KOKKOS/pair_snap_kokkos_impl.h +++ b/src/KOKKOS/pair_snap_kokkos_impl.h @@ -578,15 +578,21 @@ void PairSNAPKokkos::coeff(int narg, char d_radelem = Kokkos::View("pair:radelem",nelements); d_wjelem = Kokkos::View("pair:wjelem",nelements); d_coeffelem = Kokkos::View("pair:coeffelem",nelements,ncoeffall); + d_sinnerelem = Kokkos::View("pair:sinnerelem",nelements); + d_dinnerelem = Kokkos::View("pair:dinnerelem",nelements); auto h_radelem = Kokkos::create_mirror_view(d_radelem); auto h_wjelem = Kokkos::create_mirror_view(d_wjelem); auto h_coeffelem = Kokkos::create_mirror_view(d_coeffelem); + auto h_sinnerelem = Kokkos::create_mirror_view(d_sinnerelem); + auto h_dinnerelem = Kokkos::create_mirror_view(d_dinnerelem); auto h_map = Kokkos::create_mirror_view(d_map); for (int ielem = 0; ielem < nelements; ielem++) { h_radelem(ielem) = radelem[ielem]; h_wjelem(ielem) = wjelem[ielem]; + h_sinnerelem(ielem) = sinnerelem[ielem]; + h_dinnerelem(ielem) = dinnerelem[ielem]; for (int jcoeff = 0; jcoeff < ncoeffall; jcoeff++) { h_coeffelem(ielem,jcoeff) = coeffelem[ielem][jcoeff]; } @@ -599,10 +605,12 @@ void PairSNAPKokkos::coeff(int narg, char Kokkos::deep_copy(d_radelem,h_radelem); Kokkos::deep_copy(d_wjelem,h_wjelem); Kokkos::deep_copy(d_coeffelem,h_coeffelem); + Kokkos::deep_copy(d_sinnerelem,h_sinnerelem); + Kokkos::deep_copy(d_dinnerelem,h_dinnerelem); Kokkos::deep_copy(d_map,h_map); snaKK = SNAKokkos(rfac0,twojmax, - rmin0,switchflag,bzeroflag,chemflag,bnormflag,wselfallflag,nelements); + rmin0,switchflag,bzeroflag,chemflag,bnormflag,wselfallflag,nelements,switchinnerflag); snaKK.grow_rij(0,0); snaKK.init(); } @@ -724,15 +732,19 @@ void PairSNAPKokkos::operator() (TagPairSN const F_FLOAT dx = x(j,0) - xtmp; const F_FLOAT dy = x(j,1) - ytmp; const F_FLOAT dz = x(j,2) - ztmp; - const int elem_j = d_map[jtype]; + const int jelem = d_map[jtype]; my_sna.rij(ii,offset,0) = static_cast(dx); my_sna.rij(ii,offset,1) = static_cast(dy); my_sna.rij(ii,offset,2) = static_cast(dz); - my_sna.wj(ii,offset) = static_cast(d_wjelem[elem_j]); - my_sna.rcutij(ii,offset) = static_cast((radi + d_radelem[elem_j])*rcutfac); + my_sna.wj(ii,offset) = static_cast(d_wjelem[jelem]); + my_sna.rcutij(ii,offset) = static_cast((radi + d_radelem[jelem])*rcutfac); my_sna.inside(ii,offset) = j; + if (switchinnerflag) { + my_sna.sinnerij(ii,offset) = 0.5*(d_sinnerelem[ielem] + d_sinnerelem[jelem]); + my_sna.dinnerij(ii,offset) = 0.5*(d_dinnerelem[ielem] + d_dinnerelem[jelem]); + } if (chemflag) - my_sna.element(ii,offset) = elem_j; + my_sna.element(ii,offset) = jelem; else my_sna.element(ii,offset) = 0; } @@ -1080,18 +1092,22 @@ void PairSNAPKokkos::operator() (TagPairSN const int jtype = type(j); const F_FLOAT rsq = dx*dx + dy*dy + dz*dz; - const int elem_j = d_map[jtype]; + const int jelem = d_map[jtype]; if (rsq < rnd_cutsq(itype,jtype)) { if (final) { my_sna.rij(ii,offset,0) = static_cast(dx); my_sna.rij(ii,offset,1) = static_cast(dy); my_sna.rij(ii,offset,2) = static_cast(dz); - my_sna.wj(ii,offset) = static_cast(d_wjelem[elem_j]); - my_sna.rcutij(ii,offset) = static_cast((radi + d_radelem[elem_j])*rcutfac); + my_sna.wj(ii,offset) = static_cast(d_wjelem[jelem]); + my_sna.rcutij(ii,offset) = static_cast((radi + d_radelem[jelem])*rcutfac); my_sna.inside(ii,offset) = j; + if (switchinnerflag) { + my_sna.sinnerij(ii,offset) = 0.5*(d_sinnerelem[ielem] + d_sinnerelem[jelem]); + my_sna.dinnerij(ii,offset) = 0.5*(d_dinnerelem[ielem] + d_dinnerelem[jelem]); + } if (chemflag) - my_sna.element(ii,offset) = elem_j; + my_sna.element(ii,offset) = jelem; else my_sna.element(ii,offset) = 0; } diff --git a/src/KOKKOS/sna_kokkos.h b/src/KOKKOS/sna_kokkos.h index b3389f2a33..fe70129660 100644 --- a/src/KOKKOS/sna_kokkos.h +++ b/src/KOKKOS/sna_kokkos.h @@ -64,7 +64,7 @@ struct alignas(8) FullHalfMapper { template class SNAKokkos { -public: + public: using real_type = real_type_; using complex = SNAComplex; static constexpr int vector_length = vector_length_; @@ -97,21 +97,22 @@ public: typedef Kokkos::View t_sna_3c3; typedef Kokkos::View t_sna_5c; -inline + inline SNAKokkos() {}; + KOKKOS_INLINE_FUNCTION SNAKokkos(const SNAKokkos& sna, const typename Kokkos::TeamPolicy::member_type& team); -inline - SNAKokkos(real_type, int, real_type, int, int, int, int, int, int); + inline + SNAKokkos(real_type, int, real_type, int, int, int, int, int, int, int); KOKKOS_INLINE_FUNCTION ~SNAKokkos(); -inline + inline void build_indexlist(); // SNAKokkos() -inline + inline void init(); // double memory_usage(); @@ -192,13 +193,13 @@ inline void compute_deidrj_cpu(const typename Kokkos::TeamPolicy::member_type& team, int, int); // ForceSNAP KOKKOS_INLINE_FUNCTION - real_type compute_sfac(real_type, real_type); // add_uarraytot, compute_duarray + real_type compute_sfac(real_type, real_type, real_type, real_type); // add_uarraytot, compute_duarray KOKKOS_INLINE_FUNCTION - real_type compute_dsfac(real_type, real_type); // compute_duarray + real_type compute_dsfac(real_type, real_type, real_type, real_type); // compute_duarray KOKKOS_INLINE_FUNCTION - void compute_s_dsfac(const real_type, const real_type, real_type&, real_type&); // compute_cayley_klein + void compute_s_dsfac(const real_type, const real_type, const real_type, const real_type, real_type&, real_type&); // compute_cayley_klein #ifdef TIMING_INFO double* timers; @@ -214,6 +215,8 @@ inline t_sna_2i inside; t_sna_2d wj; t_sna_2d rcutij; + t_sna_2d sinnerij; + t_sna_2d dinnerij; t_sna_2i element; t_sna_3d dedr; int natom, nmax; @@ -255,7 +258,7 @@ inline int ndoubles; int ntriples; -private: + private: real_type rmin0, rfac0; //use indexlist instead of loops, constructor generates these @@ -264,13 +267,13 @@ private: Kokkos::View idxb; Kokkos::View idxcg_block; -public: + public: Kokkos::View idxu_block; Kokkos::View idxu_half_block; Kokkos::View idxu_cache_block; Kokkos::View idxu_full_half; -private: + private: Kokkos::View idxz_block; Kokkos::View idxb_block; @@ -290,14 +293,14 @@ private: KOKKOS_INLINE_FUNCTION void create_thread_scratch_arrays(const typename Kokkos::TeamPolicy::member_type& team); // SNAKokkos() -inline + inline void init_clebsch_gordan(); // init() -inline + inline void init_rootpqarray(); // init() KOKKOS_INLINE_FUNCTION - void add_uarraytot(const typename Kokkos::TeamPolicy::member_type& team, int, int, const real_type&, const real_type&, const real_type&, int); // compute_ui + void add_uarraytot(const typename Kokkos::TeamPolicy::member_type& team, int, int, const real_type&, const real_type&, const real_type&, const real_type&, const real_type&, int); // compute_ui KOKKOS_INLINE_FUNCTION void compute_uarray_cpu(const typename Kokkos::TeamPolicy::member_type& team, int, int, @@ -308,18 +311,24 @@ inline inline double deltacg(int, int, int); // init_clebsch_gordan -inline + inline int compute_ncoeff(); // SNAKokkos() KOKKOS_INLINE_FUNCTION void compute_duarray_cpu(const typename Kokkos::TeamPolicy::member_type& team, int, int, const real_type&, const real_type&, const real_type&, // compute_duidrj_cpu - const real_type&, const real_type&, const real_type&, const real_type&, const real_type&); + const real_type&, const real_type&, const real_type&, const real_type&, const real_type&, + const real_type&, const real_type&); // Sets the style for the switching function // 0 = none // 1 = cosine int switch_flag; + // Sets the style for the inner switching function + // 0 = none + // 1 = cosine + int switch_inner_flag; + // Chem snap flags int chem_flag; int bnorm_flag; diff --git a/src/KOKKOS/sna_kokkos_impl.h b/src/KOKKOS/sna_kokkos_impl.h index 550548df55..77130b9781 100644 --- a/src/KOKKOS/sna_kokkos_impl.h +++ b/src/KOKKOS/sna_kokkos_impl.h @@ -25,12 +25,13 @@ namespace LAMMPS_NS { static const double MY_PI = 3.14159265358979323846; // pi +static const double MY_PI2 = 1.57079632679489661923; // pi/2 template inline SNAKokkos::SNAKokkos(real_type rfac0_in, int twojmax_in, real_type rmin0_in, int switch_flag_in, int bzero_flag_in, - int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in) + int chem_flag_in, int bnorm_flag_in, int wselfall_flag_in, int nelements_in, int switch_inner_flag_in) { LAMMPS_NS::ExecutionSpace execution_space = ExecutionSpaceFromDevice::space; host_flag = (execution_space == LAMMPS_NS::Host); @@ -40,6 +41,7 @@ SNAKokkos::SNAKokkos(real_type rfac0_in, rfac0 = rfac0_in; rmin0 = rmin0_in; switch_flag = switch_flag_in; + switch_inner_flag = switch_inner_flag_in; bzero_flag = bzero_flag_in; chem_flag = chem_flag_in; @@ -295,6 +297,8 @@ void SNAKokkos::grow_rij(int newnatom, int rij = t_sna_3d(Kokkos::NoInit("sna:rij"),natom,nmax,3); wj = t_sna_2d(Kokkos::NoInit("sna:wj"),natom,nmax); rcutij = t_sna_2d(Kokkos::NoInit("sna:rcutij"),natom,nmax); + sinnerij = t_sna_2d(Kokkos::NoInit("sna:sinnerij"),natom,nmax); + dinnerij = t_sna_2d(Kokkos::NoInit("sna:dinnerij"),natom,nmax); inside = t_sna_2i(Kokkos::NoInit("sna:inside"),natom,nmax); element = t_sna_2i(Kokkos::NoInit("sna:element"),natom,nmax); dedr = t_sna_3d(Kokkos::NoInit("sna:dedr"),natom,nmax,3); @@ -371,6 +375,8 @@ void SNAKokkos::compute_cayley_klein(const const real_type rsq = x * x + y * y + z * z; const real_type r = sqrt(rsq); const real_type rcut = rcutij(iatom, jnbor); + const real_type sinner = sinnerij(iatom, jnbor); + const real_type dinner = dinnerij(iatom, jnbor); const real_type rscale0 = rfac0 * static_cast(MY_PI) / (rcut - rmin0); const real_type theta0 = (r - rmin0) * rscale0; const real_type sn = sin(theta0); @@ -380,7 +386,7 @@ void SNAKokkos::compute_cayley_klein(const const real_type wj_local = wj(iatom, jnbor); real_type sfac, dsfac; - compute_s_dsfac(r, rcut, sfac, dsfac); + compute_s_dsfac(r, rcut, sinner, dinner, sfac, dsfac); sfac *= wj_local; dsfac *= wj_local; @@ -1238,7 +1244,7 @@ void SNAKokkos::compute_ui_cpu(const typen z0 = r / tan(theta0); compute_uarray_cpu(team, iatom, jnbor, x, y, z, z0, r); - add_uarraytot(team, iatom, jnbor, r, wj(iatom,jnbor), rcutij(iatom,jnbor), element(iatom, jnbor)); + add_uarraytot(team, iatom, jnbor, r, wj(iatom,jnbor), rcutij(iatom,jnbor), sinnerij(iatom,jnbor), dinnerij(iatom,jnbor), element(iatom, jnbor)); } /* ---------------------------------------------------------------------- @@ -1541,7 +1547,7 @@ void SNAKokkos::compute_duidrj_cpu(const t z0 = r * cs / sn; dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq; - compute_duarray_cpu(team, iatom, jnbor, x, y, z, z0, r, dz0dr, wj(iatom,jnbor), rcutij(iatom,jnbor)); + compute_duarray_cpu(team, iatom, jnbor, x, y, z, z0, r, dz0dr, wj(iatom,jnbor), rcutij(iatom,jnbor), sinnerij(iatom,jnbor), dinnerij(iatom,jnbor)); } @@ -1623,9 +1629,10 @@ void SNAKokkos::compute_deidrj_cpu(const t template KOKKOS_INLINE_FUNCTION void SNAKokkos::add_uarraytot(const typename Kokkos::TeamPolicy::member_type& team, int iatom, int jnbor, - const real_type& r, const real_type& wj, const real_type& rcut, int jelem) + const real_type& r, const real_type& wj, const real_type& rcut, + const real_type& sinner, const real_type& dinner, int jelem) { - const real_type sfac = compute_sfac(r, rcut) * wj; + const real_type sfac = compute_sfac(r, rcut, sinner, dinner) * wj; Kokkos::parallel_for(Kokkos::ThreadVectorRange(team,twojmax+1), [&] (const int& j) { @@ -1746,7 +1753,8 @@ KOKKOS_INLINE_FUNCTION void SNAKokkos::compute_duarray_cpu(const typename Kokkos::TeamPolicy::member_type& team, int iatom, int jnbor, const real_type& x, const real_type& y, const real_type& z, const real_type& z0, const real_type& r, const real_type& dz0dr, - const real_type& wj, const real_type& rcut) + const real_type& wj, const real_type& rcut, + const real_type& sinner, const real_type& dinner) { real_type r0inv; real_type a_r, a_i, b_r, b_i; @@ -1872,8 +1880,8 @@ void SNAKokkos::compute_duarray_cpu(const }); } - real_type sfac = compute_sfac(r, rcut); - real_type dsfac = compute_dsfac(r, rcut); + real_type sfac = compute_sfac(r, rcut, sinner, dinner); + real_type dsfac = compute_dsfac(r, rcut, sinner, dinner); sfac *= wj; dsfac *= wj; @@ -2224,63 +2232,121 @@ int SNAKokkos::compute_ncoeff() template KOKKOS_INLINE_FUNCTION -real_type SNAKokkos::compute_sfac(real_type r, real_type rcut) +real_type SNAKokkos::compute_sfac(real_type r, real_type rcut, real_type sinner, real_type dinner) { + real_type sfac_outer; constexpr real_type one = static_cast(1.0); constexpr real_type zero = static_cast(0.0); constexpr real_type onehalf = static_cast(0.5); - if (switch_flag == 0) return one; + if (switch_flag == 0) sfac_outer = one; if (switch_flag == 1) { - if (r <= rmin0) return one; + if (r <= rmin0) sfac_outer = one; else if (r > rcut) return zero; else { real_type rcutfac = static_cast(MY_PI) / (rcut - rmin0); - return onehalf * (cos((r - rmin0) * rcutfac) + one); + sfac_outer = onehalf * (cos((r - rmin0) * rcutfac) + one); } } - return zero; + + if (switch_inner_flag == 0) return sfac_outer; + if (switch_inner_flag == 1) { + if (r >= sinner + dinner) + return sfac_outer; + else if (r > sinner - dinner) { + real_type rcutfac = static_cast(MY_PI2) / dinner; + return sfac_outer * + onehalf * (one - cos(static_cast(MY_PI2) + (r - sinner) * rcutfac)); + } else return zero; + } + return zero; // dummy return } /* ---------------------------------------------------------------------- */ template KOKKOS_INLINE_FUNCTION -real_type SNAKokkos::compute_dsfac(real_type r, real_type rcut) +real_type SNAKokkos::compute_dsfac(real_type r, real_type rcut, real_type sinner, real_type dinner) { + real_type sfac_outer, dsfac_outer, sfac_inner, dsfac_inner; + constexpr real_type one = static_cast(1.0); constexpr real_type zero = static_cast(0.0); constexpr real_type onehalf = static_cast(0.5); - if (switch_flag == 0) return zero; + if (switch_flag == 0) dsfac_outer = zero; if (switch_flag == 1) { - if (r <= rmin0) return zero; + if (r <= rmin0) dsfac_outer = zero; else if (r > rcut) return zero; else { real_type rcutfac = static_cast(MY_PI) / (rcut - rmin0); - return -onehalf * sin((r - rmin0) * rcutfac) * rcutfac; + dsfac_outer = -onehalf * sin((r - rmin0) * rcutfac) * rcutfac; } } - return zero; + + if (switch_inner_flag == 0) return dsfac_outer; + if (switch_inner_flag == 1) { + if (r >= sinner + dinner) + return dsfac_outer; + else if (r > sinner - dinner) { + + // calculate sfac_outer + + if (switch_flag == 0) sfac_outer = one; + if (switch_flag == 1) { + if (r <= rmin0) sfac_outer = one; + else if (r > rcut) sfac_outer = zero; + else { + real_type rcutfac = static_cast(MY_PI) / (rcut - rmin0); + sfac_outer = onehalf * (cos((r - rmin0) * rcutfac) + one); + } + } + + // calculate sfac_inner + + real_type rcutfac = static_cast(MY_PI2) / dinner; + sfac_inner = onehalf * (one - cos(static_cast(MY_PI2) + (r - sinner) * rcutfac)); + dsfac_inner = onehalf * rcutfac * sin(static_cast(MY_PI2) + (r - sinner) * rcutfac); + return dsfac_outer * sfac_inner + sfac_outer * dsfac_inner; + + } else return zero; + } + return zero; // dummy return } template KOKKOS_INLINE_FUNCTION -void SNAKokkos::compute_s_dsfac(const real_type r, const real_type rcut, real_type& sfac, real_type& dsfac) { +void SNAKokkos::compute_s_dsfac(const real_type r, const real_type rcut, const real_type sinner, const real_type dinner, real_type& sfac, real_type& dsfac) { + + real_type sfac_outer, dsfac_outer, sfac_inner, dsfac_inner; constexpr real_type one = static_cast(1.0); constexpr real_type zero = static_cast(0.0); constexpr real_type onehalf = static_cast(0.5); - if (switch_flag == 0) { sfac = zero; dsfac = zero; } + + if (switch_flag == 0) { sfac_outer = zero; dsfac_outer = zero; } else if (switch_flag == 1) { - if (r <= rmin0) { sfac = one; dsfac = zero; } - else if (r > rcut) { sfac = zero; dsfac = zero; } + if (r <= rmin0) { sfac_outer = one; dsfac_outer = zero; } + else if (r > rcut) { sfac = zero; dsfac = zero; return; } else { const real_type rcutfac = static_cast(MY_PI) / (rcut - rmin0); const real_type theta0 = (r - rmin0) * rcutfac; const real_type sn = sin(theta0); const real_type cs = cos(theta0); - sfac = onehalf * (cs + one); - dsfac = -onehalf * sn * rcutfac; - + sfac_outer = onehalf * (cs + one); + dsfac_outer = -onehalf * sn * rcutfac; } - } else { sfac = zero; dsfac = zero; } + } else { sfac = zero; dsfac = zero; return; } // dummy return + + if (switch_inner_flag == 0) { sfac = sfac_outer; dsfac = dsfac_outer; return; } + else if (switch_inner_flag == 1) { + if (r >= sinner + dinner) { sfac = sfac_outer; dsfac = dsfac_outer; return; } + else if (r > sinner - dinner) { + real_type rcutfac = static_cast(MY_PI2) / dinner; + sfac_inner = onehalf * (one - cos(static_cast(MY_PI2) + (r - sinner) * rcutfac)); + dsfac_inner = onehalf * rcutfac * sin(static_cast(MY_PI2) + (r - sinner) * rcutfac); + sfac = sfac_outer * sfac_inner; + dsfac = dsfac_outer * sfac_inner + sfac_outer * dsfac_inner; + return; + } else { sfac = zero; dsfac = zero; return; } + } else { sfac = zero; dsfac = zero; return; } // dummy return + } /* ---------------------------------------------------------------------- @@ -2356,6 +2422,8 @@ double SNAKokkos::memory_usage() bytes += natom * nmax * sizeof(real_type); // inside bytes += natom * nmax * sizeof(real_type); // wj bytes += natom * nmax * sizeof(real_type); // rcutij + bytes += natom * nmax * sizeof(real_type); // sinnerij + bytes += natom * nmax * sizeof(real_type); // dinnerij return bytes; } diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index b2568d3a4d..fc65ca5207 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -41,8 +41,8 @@ MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *_lmp, char *paramfilename) : ML radelem = nullptr; wjelem = nullptr; snaptr = nullptr; - rinnerelem = nullptr; - drinnerelem = nullptr; + sinnerelem = nullptr; + dinnerelem = nullptr; read_paramfile(paramfilename); @@ -59,11 +59,8 @@ MLIAPDescriptorSNAP::~MLIAPDescriptorSNAP() memory->destroy(radelem); memory->destroy(wjelem); delete snaptr; - - if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); - } + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); } /* ---------------------------------------------------------------------- @@ -94,8 +91,8 @@ void MLIAPDescriptorSNAP::compute_descriptors(class MLIAPData *data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + snaptr->sinnerij[ninside] = 0.5 * (sinnerelem[ielem] + sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5 * (dinnerelem[ielem] + dinnerelem[jelem]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -150,8 +147,8 @@ void MLIAPDescriptorSNAP::compute_forces(class MLIAPData *data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + snaptr->sinnerij[ninside] = 0.5 * (sinnerelem[ielem] + sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5 * (dinnerelem[ielem] + dinnerelem[jelem]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -222,8 +219,8 @@ void MLIAPDescriptorSNAP::compute_force_gradients(class MLIAPData *data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + snaptr->sinnerij[ninside] = 0.5 * (sinnerelem[ielem] + sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5 * (dinnerelem[ielem] + dinnerelem[jelem]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -292,8 +289,8 @@ void MLIAPDescriptorSNAP::compute_descriptor_gradients(class MLIAPData *data) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = sqrt(cutsq[ielem][jelem]); if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5 * (rinnerelem[ielem] + rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5 * (drinnerelem[ielem] + drinnerelem[jelem]); + snaptr->sinnerij[ninside] = 0.5 * (sinnerelem[ielem] + sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5 * (dinnerelem[ielem] + dinnerelem[jelem]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -364,8 +361,8 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) // set local input checks - int rinnerflag = 0; - int drinnerflag = 0; + int sinnerflag = 0; + int dinnerflag = 0; for (int i = 0; i < nelements; i++) delete[] elements[i]; delete[] elements; @@ -414,7 +411,7 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) // check for keywords with one value per element if ((keywd == "elems") || (keywd == "radelems") || (keywd == "welems") || - (keywd == "rinnerelems") || (keywd == "drinnerelems")) { + (keywd == "sinnerelems") || (keywd == "dinnerelems")) { if ((nelementsflag == 0) || ((int) words.count() != nelements + 1)) error->all(FLERR, "Incorrect SNAP parameter file"); @@ -433,14 +430,14 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) for (int ielem = 0; ielem < nelements; ielem++) wjelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); wjelemflag = 1; - } else if (keywd == "rinnerelems") { + } else if (keywd == "sinnerelems") { for (int ielem = 0; ielem < nelements; ielem++) - rinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); - rinnerflag = 1; - } else if (keywd == "drinnerelems") { + sinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); + sinnerflag = 1; + } else if (keywd == "dinnerelems") { for (int ielem = 0; ielem < nelements; ielem++) - drinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); - rinnerflag = 1; + dinnerelem[ielem] = utils::numeric(FLERR, words.next(), false, lmp); + dinnerflag = 1; } } else { @@ -457,8 +454,8 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) elements = new char *[nelements]; memory->create(radelem, nelements, "mliap_snap_descriptor:radelem"); memory->create(wjelem, nelements, "mliap_snap_descriptor:wjelem"); - memory->create(rinnerelem, nelements, "mliap_snap_descriptor:rinner"); - memory->create(drinnerelem, nelements, "mliap_snap_descriptor:drinner"); + memory->create(sinnerelem, nelements, "mliap_snap_descriptor:sinner"); + memory->create(dinnerelem, nelements, "mliap_snap_descriptor:dinner"); nelementsflag = 1; } else if (keywd == "rcutfac") { rcutfac = utils::numeric(FLERR, keyval, false, lmp); @@ -491,10 +488,10 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename) !wjelemflag) error->all(FLERR, "Incorrect SNAP parameter file"); - if (switchinnerflag && !(rinnerflag && drinnerflag)) + if (switchinnerflag && !(sinnerflag && dinnerflag)) error->all(FLERR, "Incorrect SNAP parameter file"); - if (!switchinnerflag && (rinnerflag || drinnerflag)) + if (!switchinnerflag && (sinnerflag || dinnerflag)) error->all(FLERR, "Incorrect SNAP parameter file"); // construct cutsq diff --git a/src/ML-IAP/mliap_descriptor_snap.h b/src/ML-IAP/mliap_descriptor_snap.h index 0e217f82fa..b49ff6d2cc 100644 --- a/src/ML-IAP/mliap_descriptor_snap.h +++ b/src/ML-IAP/mliap_descriptor_snap.h @@ -42,8 +42,8 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor { int switchinnerflag; double rfac0, rmin0; - double *rinnerelem; - double *drinnerelem; + double* sinnerelem; + double* dinnerelem; }; } // namespace LAMMPS_NS diff --git a/src/ML-SNAP/compute_sna_atom.cpp b/src/ML-SNAP/compute_sna_atom.cpp index cec705770d..5f24ebeaa6 100644 --- a/src/ML-SNAP/compute_sna_atom.cpp +++ b/src/ML-SNAP/compute_sna_atom.cpp @@ -32,7 +32,7 @@ using namespace LAMMPS_NS; ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), sna(nullptr), - radelem(nullptr), wjelem(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) + radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { double rmin0, rfac0; @@ -85,6 +85,11 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : } } + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + // process optional args int iarg = nargmin; @@ -134,21 +139,37 @@ ComputeSNAAtom::ComputeSNAAtom(LAMMPS *lmp, int narg, char **arg) : wselfallflag = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+1+2*ntypes > narg) + if (iarg+2 > narg) error->all(FLERR,"Illegal compute sna/atom command"); - switchinnerflag = 1; + switchinnerflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"sinner") == 0) { iarg++; - memory->create(rinnerelem,ntypes+1,"sna/atom:rinnerelem"); - memory->create(drinnerelem,ntypes+1,"sna/atom:drinnerelem"); + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute sna/atom command"); + memory->create(sinnerelem,ntypes+1,"sna/atom:sinnerelem"); for (int i = 0; i < ntypes; i++) - rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerflag = 1; iarg += ntypes; + } else if (strcmp(arg[iarg],"dinner") == 0) { + iarg++; + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute sna/atom command"); + memory->create(dinnerelem,ntypes+1,"sna/atom:dinnerelem"); for (int i = 0; i < ntypes; i++) - drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerflag = 1; iarg += ntypes; } else error->all(FLERR,"Illegal compute sna/atom command"); } + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR,"Illegal compute sna/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, @@ -176,8 +197,8 @@ ComputeSNAAtom::~ComputeSNAAtom() if (chemflag) memory->destroy(map); if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); } } @@ -282,8 +303,8 @@ void ComputeSNAAtom::compute_peratom() snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; diff --git a/src/ML-SNAP/compute_sna_atom.h b/src/ML-SNAP/compute_sna_atom.h index ffb88d4df3..7d1ebfa2f8 100644 --- a/src/ML-SNAP/compute_sna_atom.h +++ b/src/ML-SNAP/compute_sna_atom.h @@ -45,8 +45,8 @@ class ComputeSNAAtom : public Compute { int *map; // map types to [0,nelements) int nelements, chemflag; int switchinnerflag; - double *rinnerelem; - double *drinnerelem; + double *sinnerelem; + double *dinnerelem; class SNA *snaptr; double cutmax; int quadraticflag; diff --git a/src/ML-SNAP/compute_snad_atom.cpp b/src/ML-SNAP/compute_snad_atom.cpp index 4874a53c9c..838d8a85c9 100644 --- a/src/ML-SNAP/compute_snad_atom.cpp +++ b/src/ML-SNAP/compute_snad_atom.cpp @@ -32,7 +32,7 @@ using namespace LAMMPS_NS; ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snad(nullptr), - radelem(nullptr), wjelem(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) + radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; @@ -57,8 +57,8 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"sna/atom:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"sna/atom:wjelem"); + memory->create(radelem,ntypes+1,"snad/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem,ntypes+1,"snad/atom:wjelem"); rcutfac = atof(arg[3]); rfac0 = atof(arg[4]); twojmax = atoi(arg[5]); @@ -71,7 +71,7 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : double cut; cutmax = 0.0; - memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq"); + memory->create(cutsq,ntypes+1,ntypes+1,"snad/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { cut = 2.0*radelem[i]*rcutfac; if (cut > cutmax) cutmax = cut; @@ -82,6 +82,11 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : } } + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + // process optional args int iarg = nargmin; @@ -131,21 +136,38 @@ ComputeSNADAtom::ComputeSNADAtom(LAMMPS *lmp, int narg, char **arg) : wselfallflag = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+1+2*ntypes > narg) + if (iarg+2 > narg) error->all(FLERR,"Illegal compute snad/atom command"); - switchinnerflag = 1; + switchinnerflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"sinner") == 0) { iarg++; - memory->create(rinnerelem,ntypes+1,"snad/atom:rinnerelem"); - memory->create(drinnerelem,ntypes+1,"snad/atom:drinnerelem"); + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute snad/atom command"); + memory->create(sinnerelem,ntypes+1,"snad/atom:sinnerelem"); for (int i = 0; i < ntypes; i++) - rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerflag = 1; iarg += ntypes; + } else if (strcmp(arg[iarg],"dinner") == 0) { + iarg++; + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute snad/atom command"); + memory->create(dinnerelem,ntypes+1,"snad/atom:dinnerelem"); for (int i = 0; i < ntypes; i++) - drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerflag = 1; iarg += ntypes; } else error->all(FLERR,"Illegal compute snad/atom command"); } + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR,"Illegal compute snad/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, @@ -177,8 +199,8 @@ ComputeSNADAtom::~ComputeSNADAtom() if (chemflag) memory->destroy(map); if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); } } @@ -190,7 +212,7 @@ void ComputeSNADAtom::init() error->all(FLERR,"Compute snad/atom requires a pair style be defined"); if (cutmax > force->pair->cutforce) - error->all(FLERR,"Compute sna/atom cutoff is longer than pairwise cutoff"); + error->all(FLERR,"Compute snad/atom cutoff is longer than pairwise cutoff"); // need an occasional full neighbor list @@ -299,8 +321,8 @@ void ComputeSNADAtom::compute_peratom() snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; diff --git a/src/ML-SNAP/compute_snad_atom.h b/src/ML-SNAP/compute_snad_atom.h index 0951dabb0d..0cdb075daf 100644 --- a/src/ML-SNAP/compute_snad_atom.h +++ b/src/ML-SNAP/compute_snad_atom.h @@ -47,8 +47,8 @@ class ComputeSNADAtom : public Compute { int *map; // map types to [0,nelements) int nelements, chemflag; int switchinnerflag; - double *rinnerelem; - double *drinnerelem; + double *sinnerelem; + double *dinnerelem; class SNA *snaptr; double cutmax; int quadraticflag; diff --git a/src/ML-SNAP/compute_snap.cpp b/src/ML-SNAP/compute_snap.cpp index 04dac2acba..83f56e338c 100644 --- a/src/ML-SNAP/compute_snap.cpp +++ b/src/ML-SNAP/compute_snap.cpp @@ -35,7 +35,7 @@ enum{SCALAR,VECTOR,ARRAY}; ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snap(nullptr), snapall(nullptr), snap_peratom(nullptr), radelem(nullptr), wjelem(nullptr), - rinnerelem(nullptr), drinnerelem(nullptr), snaptr(nullptr) + sinnerelem(nullptr), dinnerelem(nullptr), snaptr(nullptr) { array_flag = 1; @@ -89,6 +89,11 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : } } + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + // process optional args int iarg = nargmin; @@ -115,7 +120,7 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : quadraticflag = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"chem") == 0) { - if (iarg+2 > narg) + if (iarg+2+ntypes > narg) error->all(FLERR,"Illegal compute snap command"); chemflag = 1; memory->create(map,ntypes+1,"compute_snap:map"); @@ -143,21 +148,37 @@ ComputeSnap::ComputeSnap(LAMMPS *lmp, int narg, char **arg) : bikflag = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+1+2*ntypes > narg) + if (iarg+2 > narg) error->all(FLERR,"Illegal compute snap command"); - switchinnerflag = 1; + switchinnerflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"sinner") == 0) { iarg++; - memory->create(rinnerelem,ntypes+1,"snap:rinnerelem"); - memory->create(drinnerelem,ntypes+1,"snap:drinnerelem"); + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute snap command"); + memory->create(sinnerelem,ntypes+1,"snap:sinnerelem"); for (int i = 0; i < ntypes; i++) - rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerflag = 1; iarg += ntypes; + } else if (strcmp(arg[iarg],"dinner") == 0) { + iarg++; + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute snap command"); + memory->create(dinnerelem,ntypes+1,"snap:dinnerelem"); for (int i = 0; i < ntypes; i++) - drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerflag = 1; iarg += ntypes; } else error->all(FLERR,"Illegal compute snap command"); } + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR,"Illegal compute snap command: switchinnerflag = 1, missing sinner/dinner keyword"); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR,"Illegal compute snap command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, @@ -198,8 +219,8 @@ ComputeSnap::~ComputeSnap() if (chemflag) memory->destroy(map); if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); } } @@ -353,8 +374,8 @@ void ComputeSnap::compute_array() snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; diff --git a/src/ML-SNAP/compute_snap.h b/src/ML-SNAP/compute_snap.h index 562002dcf3..bc0670e2c7 100644 --- a/src/ML-SNAP/compute_snap.h +++ b/src/ML-SNAP/compute_snap.h @@ -47,8 +47,8 @@ class ComputeSnap : public Compute { int *map; // map types to [0,nelements) int nelements, chemflag; int switchinnerflag; - double *rinnerelem; - double *drinnerelem; + double *sinnerelem; + double *dinnerelem; class SNA *snaptr; double cutmax; int quadraticflag; diff --git a/src/ML-SNAP/compute_snav_atom.cpp b/src/ML-SNAP/compute_snav_atom.cpp index d3999381f2..0baa4bc110 100644 --- a/src/ML-SNAP/compute_snav_atom.cpp +++ b/src/ML-SNAP/compute_snav_atom.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), cutsq(nullptr), list(nullptr), snav(nullptr), - radelem(nullptr), wjelem(nullptr), rinnerelem(nullptr), drinnerelem(nullptr) + radelem(nullptr), wjelem(nullptr), sinnerelem(nullptr), dinnerelem(nullptr) { double rfac0, rmin0; int twojmax, switchflag, bzeroflag, bnormflag, wselfallflag; @@ -56,8 +56,8 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : // process required arguments - memory->create(radelem,ntypes+1,"sna/atom:radelem"); // offset by 1 to match up with types - memory->create(wjelem,ntypes+1,"sna/atom:wjelem"); + memory->create(radelem,ntypes+1,"snav/atom:radelem"); // offset by 1 to match up with types + memory->create(wjelem,ntypes+1,"snav/atom:wjelem"); rcutfac = atof(arg[3]); rfac0 = atof(arg[4]); twojmax = atoi(arg[5]); @@ -67,7 +67,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : wjelem[i+1] = atof(arg[6+ntypes+i]); // construct cutsq double cut; - memory->create(cutsq,ntypes+1,ntypes+1,"sna/atom:cutsq"); + memory->create(cutsq,ntypes+1,ntypes+1,"snav/atom:cutsq"); for (int i = 1; i <= ntypes; i++) { cut = 2.0*radelem[i]*rcutfac; cutsq[i][i] = cut*cut; @@ -77,6 +77,11 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : } } + // set local input checks + + int sinnerflag = 0; + int dinnerflag = 0; + // process optional args int iarg = nargmin; @@ -104,7 +109,7 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : iarg += 2; } else if (strcmp(arg[iarg],"chem") == 0) { if (iarg+2+ntypes > narg) - error->all(FLERR,"Illegal compute sna/atom command"); + error->all(FLERR,"Illegal compute snav/atom command"); chemflag = 1; memory->create(map,ntypes+1,"compute_sna_atom:map"); nelements = utils::inumeric(FLERR,arg[iarg+1],false,lmp); @@ -126,21 +131,37 @@ ComputeSNAVAtom::ComputeSNAVAtom(LAMMPS *lmp, int narg, char **arg) : wselfallflag = atoi(arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"switchinnerflag") == 0) { - if (iarg+1+2*ntypes > narg) + if (iarg+2 > narg) error->all(FLERR,"Illegal compute snav/atom command"); - switchinnerflag = 1; + switchinnerflag = atoi(arg[iarg+1]); + iarg += 2; + } else if (strcmp(arg[iarg],"sinner") == 0) { iarg++; - memory->create(rinnerelem,ntypes+1,"snav/atom:rinnerelem"); - memory->create(drinnerelem,ntypes+1,"snav/atom:drinnerelem"); + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute snav/atom command"); + memory->create(sinnerelem,ntypes+1,"snav/atom:sinnerelem"); for (int i = 0; i < ntypes; i++) - rinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + sinnerflag = 1; iarg += ntypes; + } else if (strcmp(arg[iarg],"dinner") == 0) { + iarg++; + if (iarg+ntypes > narg) + error->all(FLERR,"Illegal compute snav/atom command"); + memory->create(dinnerelem,ntypes+1,"snav/atom:dinnerelem"); for (int i = 0; i < ntypes; i++) - drinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerelem[i+1] = utils::numeric(FLERR,arg[iarg+i],false,lmp); + dinnerflag = 1; iarg += ntypes; } else error->all(FLERR,"Illegal compute snav/atom command"); } + if (switchinnerflag && !(sinnerflag && dinnerflag)) + error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 1, missing sinner/dinner keyword"); + + if (!switchinnerflag && (sinnerflag || dinnerflag)) + error->all(FLERR,"Illegal compute snav/atom command: switchinnerflag = 0, unexpected sinner/dinner keyword"); + snaptr = new SNA(lmp, rfac0, twojmax, rmin0, switchflag, bzeroflag, chemflag, bnormflag, wselfallflag, @@ -171,8 +192,8 @@ ComputeSNAVAtom::~ComputeSNAVAtom() if (chemflag) memory->destroy(map); if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); } } @@ -291,8 +312,8 @@ void ComputeSNAVAtom::compute_peratom() snaptr->wj[ninside] = wjelem[jtype]; snaptr->rcutij[ninside] = (radi+radelem[jtype])*rcutfac; if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[itype]+rinnerelem[jtype]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[itype]+drinnerelem[jtype]); + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[itype]+sinnerelem[jtype]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[itype]+dinnerelem[jtype]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; diff --git a/src/ML-SNAP/compute_snav_atom.h b/src/ML-SNAP/compute_snav_atom.h index a261507601..1eb84d8df7 100644 --- a/src/ML-SNAP/compute_snav_atom.h +++ b/src/ML-SNAP/compute_snav_atom.h @@ -47,8 +47,8 @@ class ComputeSNAVAtom : public Compute { int *map; // map types to [0,nelements) int nelements, chemflag; int switchinnerflag; - double *rinnerelem; - double *drinnerelem; + double *sinnerelem; + double *dinnerelem; class SNA *snaptr; int quadraticflag; }; diff --git a/src/ML-SNAP/pair_snap.cpp b/src/ML-SNAP/pair_snap.cpp index 208857cd90..07a8237ab5 100644 --- a/src/ML-SNAP/pair_snap.cpp +++ b/src/ML-SNAP/pair_snap.cpp @@ -45,8 +45,8 @@ PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp) radelem = nullptr; wjelem = nullptr; coeffelem = nullptr; - rinnerelem = nullptr; - drinnerelem = nullptr; + sinnerelem = nullptr; + dinnerelem = nullptr; beta_max = 0; beta = nullptr; @@ -65,8 +65,8 @@ PairSNAP::~PairSNAP() memory->destroy(coeffelem); if (switchinnerflag) { - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); } memory->destroy(beta); @@ -158,8 +158,8 @@ void PairSNAP::compute(int eflag, int vflag) snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[ielem]+sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[ielem]+dinnerelem[jelem]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -331,8 +331,8 @@ void PairSNAP::compute_bispectrum() snaptr->wj[ninside] = wjelem[jelem]; snaptr->rcutij[ninside] = (radi + radelem[jelem])*rcutfac; if (switchinnerflag) { - snaptr->rinnerij[ninside] = 0.5*(rinnerelem[ielem]+rinnerelem[jelem]); - snaptr->drinnerij[ninside] = 0.5*(drinnerelem[ielem]+drinnerelem[jelem]); + snaptr->sinnerij[ninside] = 0.5*(sinnerelem[ielem]+sinnerelem[jelem]); + snaptr->dinnerij[ninside] = 0.5*(dinnerelem[ielem]+dinnerelem[jelem]); } if (chemflag) snaptr->element[ninside] = jelem; ninside++; @@ -516,13 +516,13 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) memory->destroy(radelem); memory->destroy(wjelem); memory->destroy(coeffelem); - memory->destroy(rinnerelem); - memory->destroy(drinnerelem); + memory->destroy(sinnerelem); + memory->destroy(dinnerelem); memory->create(radelem,nelements,"pair:radelem"); memory->create(wjelem,nelements,"pair:wjelem"); memory->create(coeffelem,nelements,ncoeffall,"pair:coeffelem"); - memory->create(rinnerelem,nelements,"pair:rinnerelem"); - memory->create(drinnerelem,nelements,"pair:drinnerelem"); + memory->create(sinnerelem,nelements,"pair:sinnerelem"); + memory->create(dinnerelem,nelements,"pair:dinnerelem"); // initialize checklist for all required nelements @@ -644,8 +644,8 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) // set local input checks - int rinnerflag = 0; - int drinnerflag = 0; + int sinnerflag = 0; + int dinnerflag = 0; // open SNAP parameter file on proc 0 @@ -688,34 +688,40 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) auto keywd = words[0]; auto keyval = words[1]; - // check for keywords with one value per element + // check for keywords with more than one value per element - if (keywd == "rinner" || keywd == "drinner") { + if (keywd == "sinner" || keywd == "dinner") { if ((int)words.size() != nelements+1) error->all(FLERR,"Incorrect SNAP parameter file"); - if (comm->me == 0) - utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, keyval); + // innerlogstr collects all values of sinner or dinner for log output below + + std::string innerlogstr; int iword = 1; - if (keywd == "rinner") { - keyval = words[iword]; + if (keywd == "sinner") { for (int ielem = 0; ielem < nelements; ielem++) { - rinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = words[iword]; + sinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); iword++; + innerlogstr += keyval + " "; } - rinnerflag = 1; - } else if (keywd == "drinner") { - keyval = words[iword]; + sinnerflag = 1; + } else if (keywd == "dinner") { for (int ielem = 0; ielem < nelements; ielem++) { - drinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); + keyval = words[iword]; + dinnerelem[ielem] = utils::numeric(FLERR,keyval,false,lmp); iword++; + innerlogstr += keyval + " "; } - drinnerflag = 1; + dinnerflag = 1; } + if (comm->me == 0) + utils::logmesg(lmp,"SNAP keyword {} {} ... \n", keywd, innerlogstr); + } else { // all other keywords take one value @@ -765,10 +771,10 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename) if (chemflag && nelemtmp != nelements) error->all(FLERR,"Incorrect SNAP parameter file"); - if (switchinnerflag && !(rinnerflag && drinnerflag)) + if (switchinnerflag && !(sinnerflag && dinnerflag)) error->all(FLERR,"Incorrect SNAP parameter file"); - if (!switchinnerflag && (rinnerflag || drinnerflag)) + if (!switchinnerflag && (sinnerflag || dinnerflag)) error->all(FLERR,"Incorrect SNAP parameter file"); } diff --git a/src/ML-SNAP/pair_snap.h b/src/ML-SNAP/pair_snap.h index 3ba492615b..fc7c41a6ce 100644 --- a/src/ML-SNAP/pair_snap.h +++ b/src/ML-SNAP/pair_snap.h @@ -60,9 +60,9 @@ class PairSNAP : public Pair { int twojmax, switchflag, bzeroflag, bnormflag; int chemflag, wselfallflag; int switchinnerflag; // inner cutoff switch - double *rinnerelem; // element inner cutoff - double *drinnerelem; // element inner cutoff range - int chunksize, parallel_thresh; + double *sinnerelem; // element inner cutoff midpoint + double *dinnerelem; // element inner cutoff half-width + int chunksize,parallel_thresh; double rfac0, rmin0, wj1, wj2; int rcutfacflag, twojmaxflag; // flags for required parameters int beta_max; // length of beta diff --git a/src/ML-SNAP/sna.cpp b/src/ML-SNAP/sna.cpp index 9d32f4ab9c..33937b9c45 100644 --- a/src/ML-SNAP/sna.cpp +++ b/src/ML-SNAP/sna.cpp @@ -145,8 +145,8 @@ SNA::SNA(LAMMPS* lmp, double rfac0_in, int twojmax_in, inside = nullptr; wj = nullptr; rcutij = nullptr; - rinnerij = nullptr; - drinnerij = nullptr; + sinnerij = nullptr; + dinnerij = nullptr; element = nullptr; nmax = 0; idxz = nullptr; @@ -176,10 +176,8 @@ SNA::~SNA() memory->destroy(inside); memory->destroy(wj); memory->destroy(rcutij); - if (switch_inner_flag) { - memory->destroy(rinnerij); - memory->destroy(drinnerij); - } + memory->destroy(sinnerij); + memory->destroy(dinnerij); if (chem_flag) memory->destroy(element); memory->destroy(ulist_r_ij); memory->destroy(ulist_i_ij); @@ -327,10 +325,8 @@ void SNA::grow_rij(int newnmax) memory->destroy(inside); memory->destroy(wj); memory->destroy(rcutij); - if (switch_inner_flag) { - memory->destroy(rinnerij); - memory->destroy(drinnerij); - } + memory->destroy(sinnerij); + memory->destroy(dinnerij); if (chem_flag) memory->destroy(element); memory->destroy(ulist_r_ij); memory->destroy(ulist_i_ij); @@ -338,10 +334,8 @@ void SNA::grow_rij(int newnmax) memory->create(inside, nmax, "pair:inside"); memory->create(wj, nmax, "pair:wj"); memory->create(rcutij, nmax, "pair:rcutij"); - if (switch_inner_flag) { - memory->create(rinnerij, nmax, "pair:rinnerij"); - memory->create(drinnerij, nmax, "pair:drinnerij"); - } + memory->create(sinnerij, nmax, "pair:sinnerij"); + memory->create(dinnerij, nmax, "pair:dinnerij"); if (chem_flag) memory->create(element, nmax, "sna:element"); memory->create(ulist_r_ij, nmax, idxu_max, "sna:ulist_ij"); memory->create(ulist_i_ij, nmax, idxu_max, "sna:ulist_ij"); @@ -1020,12 +1014,9 @@ void SNA::add_uarraytot(double r, int jj) double sfac; int jelem; - sfac = compute_sfac(r, rcutij[jj]); + sfac = compute_sfac(r, rcutij[jj], sinnerij[jj], dinnerij[jj]); sfac *= wj[jj]; - if (switch_inner_flag) - sfac *= compute_sfac_inner(r, rinnerij[jj], drinnerij[jj]); - if (chem_flag) jelem = element[jj]; else jelem = 0; @@ -1268,14 +1259,8 @@ void SNA::compute_duarray(double x, double y, double z, } } - double sfac = compute_sfac(r, rcut); - double dsfac = compute_dsfac(r, rcut); - - if (switch_inner_flag) { - double sfac_inner = compute_sfac_inner(r, rinnerij[jj], drinnerij[jj]); - dsfac = dsfac*sfac_inner + sfac*compute_dsfac_inner(r, rinnerij[jj], drinnerij[jj]); - sfac *= sfac_inner; - } + double sfac = compute_sfac(r, rcut, sinnerij[jj], dinnerij[jj]); + double dsfac = compute_dsfac(r, rcut, sinnerij[jj], dinnerij[jj]); sfac *= wj; dsfac *= wj; @@ -1339,10 +1324,8 @@ double SNA::memory_usage() bytes += (double)nmax * sizeof(int); // inside bytes += (double)nmax * sizeof(double); // wj bytes += (double)nmax * sizeof(double); // rcutij - if (switch_inner_flag) { - bytes += (double)nmax * sizeof(double); // rinnerij - bytes += (double)nmax * sizeof(double); // drinnerij - } + bytes += (double)nmax * sizeof(double); // sinnerij + bytes += (double)nmax * sizeof(double); // dinnerij if (chem_flag) bytes += (double)nmax * sizeof(int); // element return bytes; @@ -1547,9 +1530,12 @@ void SNA::compute_ncoeff() /* ---------------------------------------------------------------------- */ -double SNA::compute_sfac(double r, double rcut) +double SNA::compute_sfac(double r, double rcut, double sinner, double dinner) { double sfac; + + // calculate sfac = sfac_outer + if (switch_flag == 0) sfac = 1.0; else if (r <= rmin0) sfac = 1.0; else if (r > rcut) sfac = 0.0; @@ -1557,51 +1543,56 @@ double SNA::compute_sfac(double r, double rcut) double rcutfac = MY_PI / (rcut - rmin0); sfac = 0.5 * (cos((r - rmin0) * rcutfac) + 1.0); } + + // calculate sfac *= sfac_inner, rarely visited + + if (switch_inner_flag == 1 && r < sinner + dinner) { + if (r > sinner - dinner) { + double rcutfac = MY_PI2 / dinner; + sfac *= 0.5 * (1.0 - cos(MY_PI2 + (r - sinner) * rcutfac)); + } else sfac = 0.0; + } + return sfac; } /* ---------------------------------------------------------------------- */ -double SNA::compute_dsfac(double r, double rcut) +double SNA::compute_dsfac(double r, double rcut, double sinner, double dinner) { - double dsfac; - if (switch_flag == 0) dsfac = 0.0; - else if (r <= rmin0) dsfac = 0.0; - else if (r > rcut) dsfac = 0.0; + double dsfac, sfac_outer, dsfac_outer, sfac_inner, dsfac_inner; + if (switch_flag == 0) dsfac_outer = 0.0; + else if (r <= rmin0) dsfac_outer = 0.0; + else if (r > rcut) dsfac_outer = 0.0; else { double rcutfac = MY_PI / (rcut - rmin0); - dsfac = -0.5 * sin((r - rmin0) * rcutfac) * rcutfac; - } - return dsfac; -} - -/* ---------------------------------------------------------------------- */ - -double SNA::compute_sfac_inner(double r, double rinner, double drinner) -{ - double sfac; - if (switch_inner_flag == 0) sfac = 1.0; - else if (r >= rinner + drinner) sfac = 1.0; - else if (r <= rinner) sfac = 0.0; - else { - double rcutfac = MY_PI / drinner; - sfac = 0.5 * (1.0 - cos((r - rinner) * rcutfac)); - } - return sfac; -} - -/* ---------------------------------------------------------------------- */ - -double SNA::compute_dsfac_inner(double r, double rinner, double drinner) -{ - double dsfac; - if (switch_inner_flag == 0) dsfac = 0.0; - else if (r >= rinner + drinner) dsfac = 0.0; - else if (r <= rinner) dsfac = 0.0; - else { - double rcutfac = MY_PI / drinner; - dsfac = 0.5 * sin((r - rinner) * rcutfac) * rcutfac; + dsfac_outer = -0.5 * sin((r - rmin0) * rcutfac) * rcutfac; } + + // some duplicated computation, but rarely visited + + if (switch_inner_flag == 1 && r < sinner + dinner) { + if (r > sinner - dinner) { + + // calculate sfac_outer + + if (switch_flag == 0) sfac_outer = 1.0; + else if (r <= rmin0) sfac_outer = 1.0; + else if (r > rcut) sfac_outer = 0.0; + else { + double rcutfac = MY_PI / (rcut - rmin0); + sfac_outer = 0.5 * (cos((r - rmin0) * rcutfac) + 1.0); + } + + // calculate sfac_inner + + double rcutfac = MY_PI2 / dinner; + sfac_inner = 0.5 * (1.0 - cos(MY_PI2 + (r - sinner) * rcutfac)); + dsfac_inner = 0.5 * rcutfac * sin(MY_PI2 + (r - sinner) * rcutfac); + dsfac = dsfac_outer*sfac_inner + sfac_outer*dsfac_inner; + } else dsfac = 0.0; + } else dsfac = dsfac_outer; + return dsfac; } diff --git a/src/ML-SNAP/sna.h b/src/ML-SNAP/sna.h index b04e9dcade..d35d7bd4fd 100644 --- a/src/ML-SNAP/sna.h +++ b/src/ML-SNAP/sna.h @@ -56,11 +56,8 @@ class SNA : protected Pointers { void compute_duidrj(int); void compute_dbidrj(); void compute_deidrj(double *); - double compute_sfac(double, double); - double compute_dsfac(double, double); - - double compute_sfac_inner(double, double, double); - double compute_dsfac_inner(double, double, double); + double compute_sfac(double, double, double, double); + double compute_dsfac(double, double, double, double); // public bispectrum data @@ -80,8 +77,8 @@ class SNA : protected Pointers { // only allocated for switch_inner_flag=1 - double *rinnerij; // short inner cutoff list - double *drinnerij; // short inner range list + double *sinnerij; // short inner cutoff midpoint list + double *dinnerij; // short inner half-width list // only allocated for chem_flag=1