Merge pull request #3230 from athomps/snap-inner-mod
Snap inner cutoff improvements and porting to KOKKOS
This commit is contained in:
Axel Kohlmeyer
@ -73,7 +73,7 @@ struct TagPairSNAPComputeDeidrjCPU{};
|
||||
|
||||
template<class DeviceType, typename real_type_, int vector_length_>
|
||||
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<class DeviceType>
|
||||
inline int equal(double* x,double* y);
|
||||
inline int equal(double* x,double* y);
|
||||
template<class DeviceType>
|
||||
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<real_type*, DeviceType> d_radelem; // element radii
|
||||
Kokkos::View<real_type*, DeviceType> d_wjelem; // elements weights
|
||||
Kokkos::View<real_type**, Kokkos::LayoutRight, DeviceType> d_coeffelem; // element bispectrum coefficients
|
||||
Kokkos::View<real_type*, DeviceType> d_sinnerelem; // element inner cutoff midpoint
|
||||
Kokkos::View<real_type*, DeviceType> d_dinnerelem; // element inner cutoff half-width
|
||||
Kokkos::View<T_INT*, DeviceType> d_map; // mapping from atom types to elements
|
||||
Kokkos::View<T_INT*, DeviceType> d_ninside; // ninside for all atoms in list
|
||||
Kokkos::View<real_type**, DeviceType> d_beta; // betas for all atoms in list
|
||||
@ -329,10 +331,10 @@ inline double dist2(double* x,double* y);
|
||||
template <class DeviceType>
|
||||
class PairSNAPKokkosDevice : public PairSNAPKokkos<DeviceType, SNAP_KOKKOS_REAL, SNAP_KOKKOS_DEVICE_VECLEN> {
|
||||
|
||||
private:
|
||||
private:
|
||||
using Base = PairSNAPKokkos<DeviceType, SNAP_KOKKOS_REAL, SNAP_KOKKOS_DEVICE_VECLEN>;
|
||||
|
||||
public:
|
||||
public:
|
||||
|
||||
PairSNAPKokkosDevice(class LAMMPS *);
|
||||
|
||||
@ -348,10 +350,10 @@ public:
|
||||
template <class DeviceType>
|
||||
class PairSNAPKokkosHost : public PairSNAPKokkos<DeviceType, SNAP_KOKKOS_REAL, SNAP_KOKKOS_HOST_VECLEN> {
|
||||
|
||||
private:
|
||||
private:
|
||||
using Base = PairSNAPKokkos<DeviceType, SNAP_KOKKOS_REAL, SNAP_KOKKOS_HOST_VECLEN>;
|
||||
|
||||
public:
|
||||
public:
|
||||
|
||||
PairSNAPKokkosHost(class LAMMPS *);
|
||||
|
||||
|
||||
@ -578,15 +578,21 @@ void PairSNAPKokkos<DeviceType, real_type, vector_length>::coeff(int narg, char
|
||||
d_radelem = Kokkos::View<real_type*, DeviceType>("pair:radelem",nelements);
|
||||
d_wjelem = Kokkos::View<real_type*, DeviceType>("pair:wjelem",nelements);
|
||||
d_coeffelem = Kokkos::View<real_type**, Kokkos::LayoutRight, DeviceType>("pair:coeffelem",nelements,ncoeffall);
|
||||
d_sinnerelem = Kokkos::View<real_type*, DeviceType>("pair:sinnerelem",nelements);
|
||||
d_dinnerelem = Kokkos::View<real_type*, DeviceType>("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<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>(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<DeviceType, real_type, vector_length>::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<real_type>(dx);
|
||||
my_sna.rij(ii,offset,1) = static_cast<real_type>(dy);
|
||||
my_sna.rij(ii,offset,2) = static_cast<real_type>(dz);
|
||||
my_sna.wj(ii,offset) = static_cast<real_type>(d_wjelem[elem_j]);
|
||||
my_sna.rcutij(ii,offset) = static_cast<real_type>((radi + d_radelem[elem_j])*rcutfac);
|
||||
my_sna.wj(ii,offset) = static_cast<real_type>(d_wjelem[jelem]);
|
||||
my_sna.rcutij(ii,offset) = static_cast<real_type>((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<DeviceType, real_type, vector_length>::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<real_type>(dx);
|
||||
my_sna.rij(ii,offset,1) = static_cast<real_type>(dy);
|
||||
my_sna.rij(ii,offset,2) = static_cast<real_type>(dz);
|
||||
my_sna.wj(ii,offset) = static_cast<real_type>(d_wjelem[elem_j]);
|
||||
my_sna.rcutij(ii,offset) = static_cast<real_type>((radi + d_radelem[elem_j])*rcutfac);
|
||||
my_sna.wj(ii,offset) = static_cast<real_type>(d_wjelem[jelem]);
|
||||
my_sna.rcutij(ii,offset) = static_cast<real_type>((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;
|
||||
}
|
||||
|
||||
@ -64,7 +64,7 @@ struct alignas(8) FullHalfMapper {
|
||||
template<class DeviceType, typename real_type_, int vector_length_>
|
||||
class SNAKokkos {
|
||||
|
||||
public:
|
||||
public:
|
||||
using real_type = real_type_;
|
||||
using complex = SNAComplex<real_type>;
|
||||
static constexpr int vector_length = vector_length_;
|
||||
@ -97,21 +97,22 @@ public:
|
||||
typedef Kokkos::View<complex**[3], DeviceType> t_sna_3c3;
|
||||
typedef Kokkos::View<complex*****, DeviceType> t_sna_5c;
|
||||
|
||||
inline
|
||||
inline
|
||||
SNAKokkos() {};
|
||||
|
||||
KOKKOS_INLINE_FUNCTION
|
||||
SNAKokkos(const SNAKokkos<DeviceType,real_type,vector_length>& sna, const typename Kokkos::TeamPolicy<DeviceType>::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<DeviceType>::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<int*[3], DeviceType> idxb;
|
||||
Kokkos::View<int***, DeviceType> idxcg_block;
|
||||
|
||||
public:
|
||||
public:
|
||||
Kokkos::View<int*, DeviceType> idxu_block;
|
||||
Kokkos::View<int*, DeviceType> idxu_half_block;
|
||||
Kokkos::View<int*, DeviceType> idxu_cache_block;
|
||||
Kokkos::View<FullHalfMapper*, DeviceType> idxu_full_half;
|
||||
|
||||
private:
|
||||
private:
|
||||
Kokkos::View<int***, DeviceType> idxz_block;
|
||||
Kokkos::View<int***, DeviceType> idxb_block;
|
||||
|
||||
@ -290,14 +293,14 @@ private:
|
||||
KOKKOS_INLINE_FUNCTION
|
||||
void create_thread_scratch_arrays(const typename Kokkos::TeamPolicy<DeviceType>::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<DeviceType>::member_type& team, int, int, const real_type&, const real_type&, const real_type&, int); // compute_ui
|
||||
void add_uarraytot(const typename Kokkos::TeamPolicy<DeviceType>::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<DeviceType>::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<DeviceType>::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;
|
||||
|
||||
@ -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<class DeviceType, typename real_type, int vector_length>
|
||||
inline
|
||||
SNAKokkos<DeviceType, real_type, vector_length>::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<DeviceType>::space;
|
||||
host_flag = (execution_space == LAMMPS_NS::Host);
|
||||
@ -40,6 +41,7 @@ SNAKokkos<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>::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<real_type>(MY_PI) / (rcut - rmin0);
|
||||
const real_type theta0 = (r - rmin0) * rscale0;
|
||||
const real_type sn = sin(theta0);
|
||||
@ -380,7 +386,7 @@ void SNAKokkos<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>::compute_deidrj_cpu(const t
|
||||
template<class DeviceType, typename real_type, int vector_length>
|
||||
KOKKOS_INLINE_FUNCTION
|
||||
void SNAKokkos<DeviceType, real_type, vector_length>::add_uarraytot(const typename Kokkos::TeamPolicy<DeviceType>::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<DeviceType, real_type, vector_length>::compute_duarray_cpu(const typename Kokkos::TeamPolicy<DeviceType>::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<DeviceType, real_type, vector_length>::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<DeviceType, real_type, vector_length>::compute_ncoeff()
|
||||
|
||||
template<class DeviceType, typename real_type, int vector_length>
|
||||
KOKKOS_INLINE_FUNCTION
|
||||
real_type SNAKokkos<DeviceType, real_type, vector_length>::compute_sfac(real_type r, real_type rcut)
|
||||
real_type SNAKokkos<DeviceType, real_type, vector_length>::compute_sfac(real_type r, real_type rcut, real_type sinner, real_type dinner)
|
||||
{
|
||||
real_type sfac_outer;
|
||||
constexpr real_type one = static_cast<real_type>(1.0);
|
||||
constexpr real_type zero = static_cast<real_type>(0.0);
|
||||
constexpr real_type onehalf = static_cast<real_type>(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<real_type>(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<real_type>(MY_PI2) / dinner;
|
||||
return sfac_outer *
|
||||
onehalf * (one - cos(static_cast<real_type>(MY_PI2) + (r - sinner) * rcutfac));
|
||||
} else return zero;
|
||||
}
|
||||
return zero; // dummy return
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template<class DeviceType, typename real_type, int vector_length>
|
||||
KOKKOS_INLINE_FUNCTION
|
||||
real_type SNAKokkos<DeviceType, real_type, vector_length>::compute_dsfac(real_type r, real_type rcut)
|
||||
real_type SNAKokkos<DeviceType, real_type, vector_length>::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<real_type>(1.0);
|
||||
constexpr real_type zero = static_cast<real_type>(0.0);
|
||||
constexpr real_type onehalf = static_cast<real_type>(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<real_type>(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<real_type>(MY_PI) / (rcut - rmin0);
|
||||
sfac_outer = onehalf * (cos((r - rmin0) * rcutfac) + one);
|
||||
}
|
||||
}
|
||||
|
||||
// calculate sfac_inner
|
||||
|
||||
real_type rcutfac = static_cast<real_type>(MY_PI2) / dinner;
|
||||
sfac_inner = onehalf * (one - cos(static_cast<real_type>(MY_PI2) + (r - sinner) * rcutfac));
|
||||
dsfac_inner = onehalf * rcutfac * sin(static_cast<real_type>(MY_PI2) + (r - sinner) * rcutfac);
|
||||
return dsfac_outer * sfac_inner + sfac_outer * dsfac_inner;
|
||||
|
||||
} else return zero;
|
||||
}
|
||||
return zero; // dummy return
|
||||
}
|
||||
|
||||
template<class DeviceType, typename real_type, int vector_length>
|
||||
KOKKOS_INLINE_FUNCTION
|
||||
void SNAKokkos<DeviceType, real_type, vector_length>::compute_s_dsfac(const real_type r, const real_type rcut, real_type& sfac, real_type& dsfac) {
|
||||
void SNAKokkos<DeviceType, real_type, vector_length>::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<real_type>(1.0);
|
||||
constexpr real_type zero = static_cast<real_type>(0.0);
|
||||
constexpr real_type onehalf = static_cast<real_type>(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<real_type>(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<real_type>(MY_PI2) / dinner;
|
||||
sfac_inner = onehalf * (one - cos(static_cast<real_type>(MY_PI2) + (r - sinner) * rcutfac));
|
||||
dsfac_inner = onehalf * rcutfac * sin(static_cast<real_type>(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<DeviceType, real_type, vector_length>::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;
|
||||
}
|
||||
|
||||
@ -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
|
||||
|
||||
@ -42,8 +42,8 @@ class MLIAPDescriptorSNAP : public MLIAPDescriptor {
|
||||
int switchinnerflag;
|
||||
double rfac0, rmin0;
|
||||
|
||||
double *rinnerelem;
|
||||
double *drinnerelem;
|
||||
double* sinnerelem;
|
||||
double* dinnerelem;
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
@ -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++;
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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++;
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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++;
|
||||
|
||||
@ -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;
|
||||
|
||||
@ -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++;
|
||||
|
||||
@ -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;
|
||||
};
|
||||
|
||||
@ -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");
|
||||
}
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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;
|
||||
}
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user