diff --git a/cmake/Modules/Packages/KOKKOS.cmake b/cmake/Modules/Packages/KOKKOS.cmake index 6ac6cf651d..c3852abf7c 100644 --- a/cmake/Modules/Packages/KOKKOS.cmake +++ b/cmake/Modules/Packages/KOKKOS.cmake @@ -139,6 +139,12 @@ if(PKG_KSPACE) endif() endif() +if(PKG_ML-IAP) + list(APPEND KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/mliap_data_kokkos.cpp + ${KOKKOS_PKG_SOURCES_DIR}/mliap_descriptor_so3_kokkos.cpp + ${KOKKOS_PKG_SOURCES_DIR}/mliap_model_linear_kokkos.cpp + ${KOKKOS_PKG_SOURCES_DIR}/mliap_so3_kokkos.cpp) +endif() if(PKG_PHONON) list(APPEND KOKKOS_PKG_SOURCES ${KOKKOS_PKG_SOURCES_DIR}/dynamical_matrix_kokkos.cpp) diff --git a/doc/src/Commands_pair.rst b/doc/src/Commands_pair.rst index 90858ef9a9..23930f3f04 100644 --- a/doc/src/Commands_pair.rst +++ b/doc/src/Commands_pair.rst @@ -205,7 +205,7 @@ OPT. * :doc:`mesont/tpm ` * :doc:`mgpt ` * :doc:`mie/cut (g) ` - * :doc:`mliap ` + * :doc:`mliap (k) ` * :doc:`mm3/switch3/coulgauss/long ` * :doc:`momb ` * :doc:`morse (gkot) ` diff --git a/doc/src/pair_mliap.rst b/doc/src/pair_mliap.rst index 3993a4701e..b31ffdcea8 100644 --- a/doc/src/pair_mliap.rst +++ b/doc/src/pair_mliap.rst @@ -1,8 +1,11 @@ .. index:: pair_style mliap +.. index:: pair_style mliap/kk pair_style mliap command ======================== +Accelerator Variants: *mliap/kk* + Syntax """""" @@ -207,6 +210,12 @@ cutoff manually, such as in the following example. on the active LAMMPS object before the pair style is defined. This call locates and loads the mliap-specific python module that is built into LAMMPS. +---------- + +.. include:: accel_styles.rst + +---------- + Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index cfdf0356f7..063ea44720 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -190,6 +190,16 @@ action meam_funcs_kokkos.h meam_funcs.cpp action meam_impl_kokkos.h meam_impl.cpp action meam_setup_done_kokkos.h meam_setup_done.cpp action memory_kokkos.h +action mliap_data_kokkos.cpp mliap_data.cpp +action mliap_data_kokkos.h mliap_data.h +action mliap_descriptor_kokkos.h mliap_descriptor.h +action mliap_descriptor_so3_kokkos.cpp mliap_descriptor_so3.cpp +action mliap_descriptor_so3_kokkos.h mliap_descriptor_so3.h +action mliap_model_linear_kokkos.cpp mliap_model_linear.cpp +action mliap_model_linear_kokkos.h mliap_model_linear.h +action mliap_model_kokkos.h mliap_model.h +action mliap_so3_kokkos.cpp mliap_so3.cpp +action mliap_so3_kokkos.h mliap_so3.h action modify_kokkos.cpp action modify_kokkos.h action neigh_bond_kokkos.cpp @@ -298,6 +308,8 @@ action pair_lj_spica_kokkos.cpp pair_lj_spica.cpp action pair_lj_spica_kokkos.h pair_lj_spica.h action pair_meam_kokkos.cpp pair_meam.cpp action pair_meam_kokkos.h pair_meam.h +action pair_mliap_kokkos.cpp pair_mliap.cpp +action pair_mliap_kokkos.h pair_mliap.h action pair_morse_kokkos.cpp action pair_morse_kokkos.h action pair_multi_lucy_rx_kokkos.cpp pair_multi_lucy_rx.cpp diff --git a/src/KOKKOS/kokkos_type.h b/src/KOKKOS/kokkos_type.h index 47c50d305f..7fa07b85c5 100644 --- a/src/KOKKOS/kokkos_type.h +++ b/src/KOKKOS/kokkos_type.h @@ -709,6 +709,15 @@ typedef tdual_float_2d::t_dev_um t_float_2d_um; typedef tdual_float_2d::t_dev_const_um t_float_2d_const_um; typedef tdual_float_2d::t_dev_const_randomread t_float_2d_randomread; +//3d float array n +typedef Kokkos::DualView tdual_float_3d; +typedef tdual_float_3d::t_dev t_float_3d; +typedef tdual_float_3d::t_dev_const t_float_3d_const; +typedef tdual_float_3d::t_dev_um t_float_3d_um; +typedef tdual_float_3d::t_dev_const_um t_float_3d_const_um; +typedef tdual_float_3d::t_dev_const_randomread t_float_3d_randomread; + + //Position Types //1d X_FLOAT array n typedef Kokkos::DualView tdual_xfloat_1d; diff --git a/src/KOKKOS/memory_kokkos.h b/src/KOKKOS/memory_kokkos.h index cd163373f3..9d894a344a 100644 --- a/src/KOKKOS/memory_kokkos.h +++ b/src/KOKKOS/memory_kokkos.h @@ -175,7 +175,7 @@ TYPE create_kokkos(TYPE &data, typename TYPE::value_type **&array, array = (typename TYPE::value_type **) smalloc(nbytes,name); for (int i = 0; i < n1; i++) { - if (n2==0) + if (n2 == 0) array[i] = nullptr; else array[i] = &data.h_view(i,0); @@ -194,7 +194,7 @@ template array = (typename TYPE::value_type **) smalloc(nbytes,name); for (int i = 0; i < n1; i++) { - if (n2==0) + if (n2 == 0) array[i] = nullptr; else array[i] = &h_data(i,0); @@ -202,6 +202,67 @@ template return data; } +template + TYPE create_kokkos(TYPE &data, HTYPE &h_data, int n1, int n2, int n3, + const char *name) +{ + data = TYPE(std::string(name),n1,n2,n3); + h_data = Kokkos::create_mirror_view(data); + return data; +} + +template +TYPE create_kokkos(TYPE &data, typename TYPE::value_type ***&array, + int n1, int n2, int n3, const char *name) +{ + data = TYPE(std::string(name),n1,n2); + bigint nbytes = ((bigint) sizeof(typename TYPE::value_type **)) * n1; + array = (typename TYPE::value_type ***) smalloc(nbytes,name); + + for (int i = 0; i < n1; i++) { + if (n2 == 0) { + array[i] = nullptr; + } else { + nbytes = ((bigint) sizeof(typename TYPE::value_type *)) * n2; + array[i] = (typename TYPE::value_type **) smalloc(nbytes,name); + for (int j = 0; j < n2; j++) { + if (n3 == 0) + array[i][j] = nullptr; + else + array[i][j] = &data.h_view(i,j,0); + } + } + } + return data; +} + +template + TYPE create_kokkos(TYPE &data, HTYPE &h_data, + typename TYPE::value_type ***&array, int n1, int n2, int n3, + const char *name) +{ + data = TYPE(std::string(name),n1,n2); + h_data = Kokkos::create_mirror_view(data); + bigint nbytes = ((bigint) sizeof(typename TYPE::value_type **)) * n1; + array = (typename TYPE::value_type ***) smalloc(nbytes,name); + + for (int i = 0; i < n1; i++) { + if (n2 == 0) { + array[i] = nullptr; + } else { + nbytes = ((bigint) sizeof(typename TYPE::value_type *)) * n2; + array[i] = (typename TYPE::value_type **) smalloc(nbytes,name); + for (int j = 0; j < n2; j++) { + if (n3 == 0) + array[i][j] = nullptr; + else + array[i][j] = &data.h_view(i,j,0); + } + } + } + return data; +} + /* ---------------------------------------------------------------------- grow or shrink 1st dim of a 2d array last dim must stay the same @@ -217,7 +278,7 @@ TYPE grow_kokkos(TYPE &data, typename TYPE::value_type **&array, array = (typename TYPE::value_type**) srealloc(array,nbytes,name); for (int i = 0; i < n1; i++) - if (n2==0) + if (n2 == 0) array[i] = nullptr; else array[i] = &data.h_view(i,0); @@ -234,7 +295,7 @@ TYPE create_kokkos(TYPE &data, typename TYPE::value_type **&array, array = (typename TYPE::value_type **) smalloc(nbytes,name); for (int i = 0; i < n1; i++) - if (data.h_view.extent(1)==0) + if (data.h_view.extent(1) == 0) array[i] = nullptr; else array[i] = &data.h_view(i,0); @@ -254,7 +315,7 @@ TYPE grow_kokkos(TYPE &data, typename TYPE::value_type **&array, array = (typename TYPE::value_type **) srealloc(array,nbytes,name); for (int i = 0; i < n1; i++) - if (data.h_view.extent(1)==0) + if (data.h_view.extent(1) == 0) array[i] = nullptr; else array[i] = &data.h_view(i,0); @@ -275,6 +336,22 @@ void destroy_kokkos(TYPE data, typename TYPE::value_type** &array) array = nullptr; } +/* ---------------------------------------------------------------------- + destroy a 3d array +------------------------------------------------------------------------- */ + +template +void destroy_kokkos(TYPE data, typename TYPE::value_type*** &array) +{ + if (array == nullptr) return; + int n1 = data.extent(0); + for (int i = 0; i < n1; ++i) + sfree(array[i]); + data = TYPE(); + sfree(array); + array = nullptr; +} + /* ---------------------------------------------------------------------- reallocate Kokkos views without initialization deallocate first to reduce memory use diff --git a/src/KOKKOS/mliap_data_kokkos.cpp b/src/KOKKOS/mliap_data_kokkos.cpp new file mode 100644 index 0000000000..f9453301c7 --- /dev/null +++ b/src/KOKKOS/mliap_data_kokkos.cpp @@ -0,0 +1,341 @@ +// clang-format off +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) +------------------------------------------------------------------------- */ + +#include "mliap_data_kokkos.h" + +#include "atom_kokkos.h" +#include "kokkos_type.h" +#include "pair_mliap_kokkos.h" +#include "atom_masks.h" +#include "mliap_descriptor.h" +#include "lammps.h" +#include "kokkos.h" + +/* ---------------------------------------------------------------------- */ + +namespace LAMMPS_NS { +template +MLIAPDataKokkos::MLIAPDataKokkos(LAMMPS *lmp_in, int gradgradflag_in, int *map_in, + class MLIAPModel* model_in, + class MLIAPDescriptor* descriptor_in, + class PairMLIAPKokkos* pairmliap_in) : + MLIAPData(lmp_in, gradgradflag_in, map_in, model_in, descriptor_in, pairmliap_in), + k_pairmliap(pairmliap_in), + lmp(lmp_in) +{ + execution_space = ExecutionSpaceFromDevice::space; +} + +/* ---------------------------------------------------------------------- */ + +template +MLIAPDataKokkos::~MLIAPDataKokkos() { + memoryKK->destroy_kokkos(k_gradforce,gradforce); + memoryKK->destroy_kokkos(k_betas,betas); + memoryKK->destroy_kokkos(k_descriptors,descriptors); + memoryKK->destroy_kokkos(k_eatoms,eatoms); + memoryKK->destroy_kokkos(k_gamma_row_index,gamma_row_index); + memoryKK->destroy_kokkos(k_gamma_col_index,gamma_col_index); + memoryKK->destroy_kokkos(k_gamma,gamma); + memoryKK->destroy_kokkos(k_iatoms,iatoms); + memoryKK->destroy_kokkos(k_ielems,ielems); + memoryKK->destroy_kokkos(k_numneighs,numneighs); + memoryKK->destroy_kokkos(k_jatoms,jatoms); + memoryKK->destroy_kokkos(k_jelems,jelems); + memoryKK->destroy_kokkos(k_ij); + memoryKK->destroy_kokkos(k_rij,rij); + memoryKK->destroy_kokkos(k_graddesc,graddesc); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDataKokkos::generate_neighdata(class NeighList *list_in, int eflag_in, int vflag_in) { + + list = list_in; + + // grow nmax gradforce array if necessary + + if (atom->nmax > nmax) { + nmax = atom->nmax; + memoryKK->destroy_kokkos(k_gradforce,gradforce); + memoryKK->create_kokkos(k_gradforce, gradforce, nmax, size_gradforce, "mliap_data:gradforce"); + } + + // clear gradforce array + + auto d_gradforce = k_gradforce.template view(); + Kokkos::deep_copy(d_gradforce, 0.); + + // grow arrays if necessary + + nlistatoms = list->inum; + if (nlistatoms_max < nlistatoms) { + memoryKK->destroy_kokkos(k_betas,betas); + memoryKK->create_kokkos(k_betas, betas, nlistatoms, ndescriptors, "mliap_data:betas"); + memoryKK->destroy_kokkos(k_descriptors,descriptors); + memoryKK->create_kokkos(k_descriptors, descriptors, nlistatoms, ndescriptors, "mliap_data:descriptors"); + memoryKK->destroy_kokkos(k_eatoms,eatoms); + memoryKK->create_kokkos(k_eatoms, eatoms, nlistatoms, "mliap_data:eatoms"); + nlistatoms_max = nlistatoms; + } + + // grow gamma arrays if necessary + + if (gradgradflag == 1) { + if (natomgamma_max < nlistatoms) { + memoryKK->destroy_kokkos(k_gamma_row_index,gamma_row_index); + memoryKK->create_kokkos(k_gamma_row_index, gamma_row_index, nlistatoms, gamma_nnz, "mliap_data:gamma_row_index"); + memoryKK->destroy_kokkos(k_gamma_col_index,gamma_col_index); + memoryKK->create_kokkos(k_gamma_col_index, gamma_col_index, nlistatoms, gamma_nnz, "mliap_data:gamma_col_index"); + memoryKK->destroy_kokkos(k_gamma,gamma); + memoryKK->create_kokkos(k_gamma, gamma, nlistatoms, gamma_nnz, "mliap_data:"); + natomgamma_max = nlistatoms; + } + } + atomKK->sync(execution_space,X_MASK | V_MASK | F_MASK | TYPE_MASK); + k_pairmliap->x = atomKK->k_x.view(); + k_pairmliap->v = atomKK->k_v.view(); + k_pairmliap->f = atomKK->k_f.view(); + + grow_neigharrays(); + + // Use the ielems memory for prefix scan and set it at the end to the i type + + auto d_iatoms = k_iatoms.template view(); + auto d_ielems = k_ielems.template view(); + auto d_ij = k_ij.template view(); + auto d_numneighs = k_numneighs.template view(); + auto d_jatoms = k_jatoms.template view(); + auto d_jelems= k_jelems.template view(); + auto d_rij= k_rij.template view(); + + NeighListKokkos* k_list = static_cast*>(list); + auto d_numneigh = k_list->d_numneigh; + auto d_neighbors = k_list->d_neighbors; + auto d_ilist = k_list->d_ilist; + auto d_cutsq=k_pairmliap->k_cutsq.template view(); + + AtomKokkos *atomKK = (AtomKokkos *) atom; + auto x = atomKK->k_x.view(); + auto type = atomKK->k_type.view(); + auto map=k_pairmliap->k_map.template view(); + + Kokkos::parallel_scan(nlistatoms, KOKKOS_LAMBDA (int ii, int &update, const bool final) { + if (final) + d_ij(ii) = update; + update += d_numneighs(ii); + }); + + Kokkos::parallel_for(nlistatoms, KOKKOS_LAMBDA (int ii) { + int ij = d_ij(ii); + const int i = d_ilist[ii]; + const double xtmp = x(i, 0); + const double ytmp = x(i, 1); + const double ztmp = x(i, 2); + const int itype = type(i); + const int ielem = map(itype); + const int jnum = d_numneigh(i); + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + const double delx = x(j,0) - xtmp; + const double dely = x(j,1) - ytmp; + const double delz = x(j,2) - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type(j); + const int jelem = map(jtype); + if (rsq < d_cutsq(itype,jtype)) { + d_jatoms(ij) = j; + d_jelems(ij) = jelem; + d_rij(ij, 0) = delx; + d_rij(ij, 1) = dely; + d_rij(ij, 2) = delz; + ij++; + } + } + d_iatoms[ii] = i; + d_ielems[ii] = ielem; + }); + + modified(execution_space, NUMNEIGHS_MASK | IATOMS_MASK | IELEMS_MASK | JATOMS_MASK | JELEMS_MASK | RIJ_MASK | IJ_MASK ); + eflag = eflag_in; + vflag = vflag_in; +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDataKokkos::grow_neigharrays() { + AtomKokkos *atomKK = (AtomKokkos *) atom; + + // grow neighbor arrays if necessary + + if (natomneigh_max < nlistatoms) { + natomneigh_max = nlistatoms; + + memoryKK->destroy_kokkos(k_iatoms,iatoms); + memoryKK->create_kokkos(k_iatoms, iatoms, natomneigh_max, "mliap_data:iatoms"); + memoryKK->destroy_kokkos(k_ielems,ielems); + memoryKK->create_kokkos(k_ielems, ielems, natomneigh_max, "mliap_data:ielems"); + memoryKK->destroy_kokkos(k_ij); + memoryKK->create_kokkos(k_ij, natomneigh_max, "mliap_data:ij"); + memoryKK->destroy_kokkos(k_numneighs,numneighs); + memoryKK->create_kokkos(k_numneighs, numneighs, natomneigh_max, "mliap_data:numneighs"); + } + + NeighListKokkos* k_list = static_cast*>(list); + auto d_numneigh = k_list->d_numneigh; + auto d_neighbors = k_list->d_neighbors; + auto d_ilist = k_list->d_ilist; + + auto x = atomKK->k_x.view(); + auto type = atomKK->k_type.view(); + auto d_cutsq=k_pairmliap->k_cutsq.template view(); + auto d_numneighs = k_numneighs.template view(); + Kokkos::parallel_reduce(nlistatoms, KOKKOS_LAMBDA (int ii, int &contrib) { + const int i = d_ilist[ii]; + int count=0; + const double xtmp = x(i, 0); + const double ytmp = x(i, 1); + const double ztmp = x(i, 2); + const int itype = type(i); + const int jnum = d_numneigh(i); + for (int jj = 0; jj < jnum; jj++) { + int j = d_neighbors(i,jj); + j &= NEIGHMASK; + const double delx = x(j,0) - xtmp; + const double dely = x(j,1) - ytmp; + const double delz = x(j,2) - ztmp; + const double rsq = delx*delx + dely*dely + delz*delz; + int jtype = type(j); + if (rsq < d_cutsq(itype,jtype)) + count++; + } + d_numneighs(ii) = count; + contrib += count; + }, nij_total); + modified(execution_space, NUMNEIGHS_MASK); + + if (nneigh_max < nij_total) { + memoryKK->destroy_kokkos(k_jatoms,jatoms); + memoryKK->create_kokkos(k_jatoms, jatoms, nij_total, "mliap_data:jatoms"); + memoryKK->destroy_kokkos(k_jelems,jelems); + memoryKK->create_kokkos(k_jelems, jelems, nij_total, "mliap_data:jelems"); + memoryKK->destroy_kokkos(k_rij,rij); + memoryKK->create_kokkos(k_rij, rij, nij_total, 3, "mliap_data:rij"); + + if (gradgradflag == 0){ + memoryKK->destroy_kokkos(k_graddesc,graddesc); + memoryKK->create_kokkos(k_graddesc, graddesc, nij_total, ndescriptors,3, "mliap_data:graddesc"); + } + nneigh_max = nij_total; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDataKokkos::modified(ExecutionSpace space, unsigned int mask, bool ignore_auto_sync) { + if (space == Device) { + if (mask & IATOMS_MASK ) k_iatoms .modify(); + if (mask & IELEMS_MASK ) k_ielems .modify(); + if (mask & JATOMS_MASK ) k_jatoms .modify(); + if (mask & JELEMS_MASK ) k_jelems .modify(); + if (mask & IJ_MASK ) k_ij .modify(); + if (mask & BETAS_MASK ) k_betas .modify(); + if (mask & DESCRIPTORS_MASK ) k_descriptors .modify(); + if (mask & EATOMS_MASK ) k_eatoms .modify(); + if (mask & RIJ_MASK ) k_rij .modify(); + if (mask & GRADFORCE_MASK ) k_gradforce .modify(); + if (mask & GRADDESC_MASK ) k_graddesc .modify(); + if (mask & NUMNEIGHS_MASK ) k_numneighs .modify(); + if (mask & GAMMA_MASK_MASK ) k_gamma .modify(); + if (mask & GAMMA_ROW_MASK ) k_gamma_row_index.modify(); + if (mask & GAMMA_COL_MASK ) k_gamma_col_index.modify(); + + if (lmp->kokkos->auto_sync && !ignore_auto_sync) sync(Host, mask, true); + } else { + if (mask & IATOMS_MASK ) k_iatoms .modify(); + if (mask & IELEMS_MASK ) k_ielems .modify(); + if (mask & JATOMS_MASK ) k_jatoms .modify(); + if (mask & JELEMS_MASK ) k_jelems .modify(); + if (mask & IJ_MASK ) k_ij .modify(); + if (mask & BETAS_MASK ) k_betas .modify(); + if (mask & DESCRIPTORS_MASK ) k_descriptors .modify(); + if (mask & EATOMS_MASK ) k_eatoms .modify(); + if (mask & RIJ_MASK ) k_rij .modify(); + if (mask & GRADFORCE_MASK ) k_gradforce .modify(); + if (mask & GRADDESC_MASK ) k_graddesc .modify(); + if (mask & NUMNEIGHS_MASK ) k_numneighs .modify(); + if (mask & GAMMA_MASK_MASK ) k_gamma .modify(); + if (mask & GAMMA_ROW_MASK ) k_gamma_row_index.modify(); + if (mask & GAMMA_COL_MASK ) k_gamma_col_index.modify(); + if (lmp->kokkos->auto_sync && !ignore_auto_sync) sync(Device, mask, true); + } +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDataKokkos::sync(ExecutionSpace space, unsigned int mask, bool ignore_auto_sync) { + + if (space == Device) { + if (lmp->kokkos->auto_sync && !ignore_auto_sync) modified(Host, mask, true); + if (mask & IATOMS_MASK ) k_iatoms .sync(); + if (mask & IELEMS_MASK ) k_ielems .sync(); + if (mask & JATOMS_MASK ) k_jatoms .sync(); + if (mask & JELEMS_MASK ) k_jelems .sync(); + if (mask & IJ_MASK ) k_ij .sync(); + if (mask & BETAS_MASK ) k_betas .sync(); + if (mask & DESCRIPTORS_MASK ) k_descriptors .sync(); + if (mask & EATOMS_MASK ) k_eatoms .sync(); + if (mask & RIJ_MASK ) k_rij .sync(); + if (mask & GRADFORCE_MASK ) k_gradforce .sync(); + if (mask & GRADDESC_MASK ) k_graddesc .sync(); + if (mask & NUMNEIGHS_MASK ) k_numneighs .sync(); + if (mask & GAMMA_MASK_MASK ) k_gamma .sync(); + if (mask & GAMMA_ROW_MASK ) k_gamma_row_index.sync(); + if (mask & GAMMA_COL_MASK ) k_gamma_col_index.sync(); + } else { + if (lmp->kokkos->auto_sync && !ignore_auto_sync) modified(Device, mask, true); + if (mask & IATOMS_MASK ) k_iatoms .sync(); + if (mask & IELEMS_MASK ) k_ielems .sync(); + if (mask & JATOMS_MASK ) k_jatoms .sync(); + if (mask & JELEMS_MASK ) k_jelems .sync(); + if (mask & IJ_MASK ) k_ij .sync(); + if (mask & BETAS_MASK ) k_betas .sync(); + if (mask & DESCRIPTORS_MASK ) k_descriptors .sync(); + if (mask & EATOMS_MASK ) k_eatoms .sync(); + if (mask & RIJ_MASK ) k_rij .sync(); + if (mask & GRADFORCE_MASK ) k_gradforce .sync(); + if (mask & GRADDESC_MASK ) k_graddesc .sync(); + if (mask & NUMNEIGHS_MASK ) k_numneighs .sync(); + if (mask & GAMMA_MASK_MASK ) k_gamma .sync(); + if (mask & GAMMA_ROW_MASK ) k_gamma_row_index.sync(); + if (mask & GAMMA_COL_MASK ) k_gamma_col_index.sync(); + } +} + +/* ---------------------------------------------------------------------- */ + +template class MLIAPDataKokkos; +#ifdef LMP_KOKKOS_GPU +template class MLIAPDataKokkos; +#endif +}// namespace diff --git a/src/KOKKOS/mliap_data_kokkos.h b/src/KOKKOS/mliap_data_kokkos.h new file mode 100644 index 0000000000..ba81e2a226 --- /dev/null +++ b/src/KOKKOS/mliap_data_kokkos.h @@ -0,0 +1,87 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) +------------------------------------------------------------------------- */ + +#ifndef LMP_MLIAP_DATA_KOKKOS_H +#define LMP_MLIAP_DATA_KOKKOS_H + +#include "mliap_data.h" + +#include "kokkos_type.h" +#include "memory_kokkos.h" +#include "pair_mliap_kokkos.h" +#include "pointers.h" + +namespace LAMMPS_NS { +// clang-format off +enum { + IATOMS_MASK = 0x00000001, + IELEMS_MASK = 0x00000002, + JATOMS_MASK = 0x00000004, + JELEMS_MASK = 0x00000008, + IJ_MASK = 0x00000010, + BETAS_MASK = 0x00000020, + DESCRIPTORS_MASK = 0x00000040, + EATOMS_MASK = 0x00000080, + RIJ_MASK = 0x00000100, + GRADFORCE_MASK = 0x00000200, + GRADDESC_MASK = 0x00000400, + NUMNEIGHS_MASK = 0x00000800, + GAMMA_MASK_MASK = 0x00001000, + GAMMA_ROW_MASK = 0x00002000, + GAMMA_COL_MASK = 0x00004000, +}; +// clang-format on + +template class MLIAPDataKokkos : public MLIAPData { + public: + MLIAPDataKokkos(class LAMMPS *, int, int *, class MLIAPModel *, class MLIAPDescriptor *, + class PairMLIAPKokkos * = nullptr); + ~MLIAPDataKokkos() override; + ExecutionSpace execution_space; + + void generate_neighdata(class NeighList *, int = 0, int = 0) override; + void grow_neigharrays() override; + + void modified(ExecutionSpace space, unsigned int mask, bool ignore_auto_sync = false); + + void sync(ExecutionSpace space, unsigned int mask, bool ignore_auto_sync = false); + + PairMLIAPKokkos *k_pairmliap; + + DAT::tdual_int_1d k_iatoms; // index of each atom + DAT::tdual_int_1d k_ielems; // element of each atom + DAT::tdual_int_1d k_jatoms; // index of each neighbor + DAT::tdual_int_1d k_jelems; // element of each neighbor + DAT::tdual_int_1d k_ij; // Start location for each particle + DAT::tdual_float_2d k_betas; // betas for all atoms in list + DAT::tdual_float_2d k_descriptors; // descriptors for all atoms in list + DAT::tdual_float_1d k_eatoms; // energies for all atoms in list + DAT::tdual_float_2d k_rij; // distance vector of each neighbor + DAT::tdual_float_2d k_gradforce; + DAT::tdual_float_3d k_graddesc; // descriptor gradient w.r.t. each neighbor + DAT::tdual_int_1d k_numneighs; // neighbors count for each atom + DAT::tdual_float_2d k_gamma; // gamma element + DAT::tdual_int_2d k_gamma_row_index; // row (parameter) index + DAT::tdual_int_2d k_gamma_col_index; // column (descriptor) index + + int nij_total; + + protected: + class LAMMPS *lmp; +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/KOKKOS/mliap_descriptor_kokkos.h b/src/KOKKOS/mliap_descriptor_kokkos.h new file mode 100644 index 0000000000..cb02a81648 --- /dev/null +++ b/src/KOKKOS/mliap_descriptor_kokkos.h @@ -0,0 +1,54 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) +------------------------------------------------------------------------- */ + +#ifndef LMP_MLIAP_DESCRIPTOR_KOKKOS_H +#define LMP_MLIAP_DESCRIPTOR_KOKKOS_H + +#include "kokkos_type.h" +#include "memory_kokkos.h" +#include "mliap_descriptor.h" +#include "pointers.h" + +namespace LAMMPS_NS { +template class MLIAPDescriptorKokkos : virtual protected Pointers { + public: + MLIAPDescriptorKokkos(LAMMPS *lmp, MLIAPDescriptor *descriptor_in) : + Pointers(lmp), descriptor(descriptor_in) + { + memoryKK->destroy_kokkos(k_wjelem); + } + + void init_data() + { + int num_elems = descriptor->nelements; + memoryKK->destroy_kokkos(k_wjelem); + memoryKK->create_kokkos(k_wjelem, num_elems, "MLIAPDescriptorKokkos::k_wjelem"); + for (int i = 0; i < num_elems; ++i) k_wjelem.h_view(i) = descriptor->wjelem[i]; + k_wjelem.modify(); + k_wjelem.sync(); + } + + virtual ~MLIAPDescriptorKokkos() + { + memoryKK->destroy_kokkos(k_wjelem); + } + + MLIAPDescriptor *descriptor; + DAT::tdual_float_1d k_wjelem; +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp b/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp new file mode 100644 index 0000000000..177ddff80e --- /dev/null +++ b/src/KOKKOS/mliap_descriptor_so3_kokkos.cpp @@ -0,0 +1,275 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#include "mliap_descriptor_so3_kokkos.h" + +#include "atom_kokkos.h" +#include "comm.h" +#include "error.h" +#include "memory.h" +#include "mliap_data_kokkos.h" +#include "mliap_so3_kokkos.h" +#include "pair_mliap.h" +#include "tokenizer.h" + +#include + +using namespace LAMMPS_NS; + +#define MAXLINE 1024 +#define MAXWORD 3 + +/* ---------------------------------------------------------------------- */ +template +MLIAPDescriptorSO3Kokkos::MLIAPDescriptorSO3Kokkos(LAMMPS *lmp, char *paramfilename) + // TODO: why take self as param, shouldn't be needed + : Pointers(lmp), MLIAPDescriptorSO3(lmp, paramfilename), MLIAPDescriptorKokkos(lmp, this) +{ + // TODO: the MLIAP_SO3 object likely needs a kokkos-ified version + so3ptr_kokkos = new MLIAP_SO3Kokkos(lmp, rcutfac, lmax, nmax, alpha); +} + +/* ---------------------------------------------------------------------- */ + +template +MLIAPDescriptorSO3Kokkos::~MLIAPDescriptorSO3Kokkos() +{ +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDescriptorSO3Kokkos::compute_descriptors(class MLIAPData *data_) +{ + auto data = static_cast*>(data_); + so3ptr_kokkos->spectrum(data->nlistatoms, data->k_numneighs, data->k_jelems, this->k_wjelem, data->k_rij, data->k_ij, + nmax, lmax, rcutfac, alpha, data->nij_total, data->ndescriptors); + + Kokkos::deep_copy(data->k_descriptors.template view(), so3ptr_kokkos->m_plist_r); + Kokkos::deep_copy(data->k_descriptors.h_view, so3ptr_kokkos->m_plist_r); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDescriptorSO3Kokkos::compute_forces(class MLIAPData *data_) +{ + auto data = static_cast*>(data_); + int npairs = data->nij_total; + auto d_numneighs = data->k_numneighs.template view(); + so3ptr_kokkos->spectrum_dxdr(data->nlistatoms, data->k_numneighs, data->k_jelems, this->k_wjelem, data->k_rij, data->k_ij, + nmax, lmax, rcutfac, alpha, npairs, data->ndescriptors); + + auto d_f = atomKK->k_f.view(); + auto d_iatoms = data->k_iatoms.template view(); + auto d_jatoms = data->k_jatoms.template view(); + auto d_betas = data->k_betas.template view(); + auto d_rij = data->k_rij.template view(); + auto d_ij = data->k_ij.template view(); + auto ndescriptors = data->ndescriptors; + auto d_dplist_r = so3ptr_kokkos->k_dplist_r; + auto vflag=data->vflag; + int vflag_either=data->k_pairmliap->vflag_either, vflag_global=data->pairmliap->vflag_global, vflag_atom=data->pairmliap->vflag_atom; + auto d_vatom = data->k_pairmliap->k_vatom.template view(); + Kokkos::View virial("virial"); + data->k_pairmliap->k_vatom.template modify(); + data->k_pairmliap->k_vatom.template sync(); + Kokkos::parallel_for(data->nlistatoms, KOKKOS_LAMBDA(int ii) { + double fij[3]; + const int i = d_iatoms(ii); + + // ensure rij, inside, wj, and rcutij are of size jnum + + const int jnum = d_numneighs(ii); + int ij = d_ij(ii); // use precomputed ij to break loop dependency + for (int jj = 0; jj < jnum; jj++) { + int j = d_jatoms(ij); + + for (int ir = 0; ir < 3; ir++) { + fij[ir] = 0.0; + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) { + fij[ir] += d_betas(ii, icoeff) * + d_dplist_r(ij,icoeff, ir); + } + } + + Kokkos::atomic_add(&d_f(i, 0),fij[0]); + Kokkos::atomic_add(&d_f(i, 1),fij[1]); + Kokkos::atomic_add(&d_f(i, 2),fij[2]); + Kokkos::atomic_add(&d_f(j, 0),-fij[0]); + Kokkos::atomic_add(&d_f(j, 1),-fij[1]); + Kokkos::atomic_add(&d_f(j, 2),-fij[2]); + + if (vflag) { + v_tally(vflag_either, vflag_global, vflag_atom, i, j, ij, fij, d_rij, virial, d_vatom); + } + ij++; + } + }); + + if (vflag) { + if (vflag_global) { + Kokkos::View h_virial("h_virial"); + Kokkos::deep_copy(h_virial,virial); + for (int i=0;i<6;++i) + data->k_pairmliap->virial[i]+=h_virial[i]; + } + if (vflag_atom) { + data->k_pairmliap->k_vatom.template modify(); + data->k_pairmliap->k_vatom.template sync(); + } + } +} + +/* ---------------------------------------------------------------------- + add virial contribution into global and per-atom accumulators +------------------------------------------------------------------------- */ +template +template +KOKKOS_INLINE_FUNCTION +void MLIAPDescriptorSO3Kokkos::v_tally(int vflag_either, int vflag_global, int vflag_atom, int i, int j, int ij, + double *fij, ViewType rij, Kokkos::View virial, ViewType vatom) +{ + double v[6]; + if (vflag_either) { + v[0] = -rij(ij,0)*fij[0]; + v[1] = -rij(ij,1)*fij[1]; + v[2] = -rij(ij,2)*fij[2]; + v[3] = -rij(ij,0)*fij[1]; + v[4] = -rij(ij,0)*fij[2]; + v[5] = -rij(ij,1)*fij[2]; + if (vflag_global) { + Kokkos::atomic_add(&virial[0], v[0]); + Kokkos::atomic_add(&virial[1], v[1]); + Kokkos::atomic_add(&virial[2], v[2]); + Kokkos::atomic_add(&virial[3], v[3]); + Kokkos::atomic_add(&virial[4], v[4]); + Kokkos::atomic_add(&virial[5], v[5]); + } + if (vflag_atom) { + Kokkos::atomic_add(&vatom(i,0), 0.5*v[0]); + Kokkos::atomic_add(&vatom(i,1), 0.5*v[1]); + Kokkos::atomic_add(&vatom(i,2), 0.5*v[2]); + Kokkos::atomic_add(&vatom(i,3), 0.5*v[3]); + Kokkos::atomic_add(&vatom(i,4), 0.5*v[4]); + Kokkos::atomic_add(&vatom(i,5), 0.5*v[5]); + + Kokkos::atomic_add(&vatom(j,0), 0.5*v[0]); + Kokkos::atomic_add(&vatom(j,1), 0.5*v[1]); + Kokkos::atomic_add(&vatom(j,2), 0.5*v[2]); + Kokkos::atomic_add(&vatom(j,3), 0.5*v[3]); + Kokkos::atomic_add(&vatom(j,4), 0.5*v[4]); + Kokkos::atomic_add(&vatom(j,5), 0.5*v[5]); + } + } +} +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDescriptorSO3Kokkos::compute_force_gradients(class MLIAPData *data_) +{ + error->all(FLERR,"This has not been tested in cuda/kokkos"); + + auto data = static_cast*>(data_); + int npairs = data->nij_total; + so3ptr_kokkos->spectrum_dxdr(data->nlistatoms, data->k_numneighs, data->k_jelems, this->k_wjelem, data->k_rij, data->k_ij, + nmax, lmax, rcutfac, alpha, npairs, data->ndescriptors); + auto d_dplist_r = so3ptr_kokkos->k_dplist_r; + auto d_gradforce = data->k_gradforce.template view(); + auto d_gamma = data->k_gamma.template view(); + auto d_gamma_row_index = data->k_gamma_row_index.template view(); + auto d_gamma_col_index = data->k_gamma_col_index.template view(); + auto d_jatoms = data->k_jatoms.template view(); + auto d_ij = data->k_ij.template view(); + auto d_numneighs = data->k_numneighs.template view(); + auto d_iatoms = data->k_iatoms.template view(); + + auto yoffset = data->yoffset, zoffset = data->zoffset, gamma_nnz = data->gamma_nnz; + + Kokkos::parallel_for (data->nlistatoms, KOKKOS_LAMBDA (int ii) { + const int i = d_iatoms(ii); + + // insure rij, inside, wj, and rcutij are of size jnum + + const int jnum = d_numneighs(ii); + int ij = d_ij(ii); + for (int jj = 0; jj < jnum; jj++) { + int j = d_jatoms(ij); + + for (int inz = 0; inz < gamma_nnz; inz++) { + const int l = d_gamma_row_index(ii, inz); + const int k = d_gamma_col_index(ii, inz); + + Kokkos::atomic_add(&d_gradforce(i, l), + + d_gamma(ii, inz) * d_dplist_r(ij, k, 0)); + Kokkos::atomic_add(&d_gradforce(i, l + yoffset), + + d_gamma(ii, inz) * d_dplist_r(ij, k, 1)); + Kokkos::atomic_add(&d_gradforce(i, l + zoffset), + + d_gamma(ii, inz) * d_dplist_r(ij, k, 2)); + Kokkos::atomic_add(&d_gradforce(j, l), - + d_gamma(ii, inz) * d_dplist_r(ij, k, 0)); + Kokkos::atomic_add(&d_gradforce(j, l + yoffset), - + d_gamma(ii, inz) * d_dplist_r(ij, k, 1)); + Kokkos::atomic_add(&d_gradforce(j, l + zoffset), - + d_gamma(ii, inz) * d_dplist_r(ij, k, 2)); + } + ij++; + } + }); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDescriptorSO3Kokkos::compute_descriptor_gradients(class MLIAPData *data_) +{ + auto data = static_cast*>(data_); + bigint npairs = data->nij_total; + so3ptr_kokkos->spectrum_dxdr(data->nlistatoms, data->k_numneighs, data->k_jelems, this->k_wjelem, data->k_rij, data->k_ij, + nmax, lmax, rcutfac, alpha, npairs, data->ndescriptors); + auto graddesc = data->k_graddesc.template view(); + Kokkos::deep_copy(graddesc, so3ptr_kokkos->k_dplist_r); + Kokkos::deep_copy(data->k_graddesc.h_view, graddesc); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPDescriptorSO3Kokkos::init() +{ + so3ptr_kokkos->init(); + MLIAPDescriptorKokkos::init_data(); +} + +/* ---------------------------------------------------------------------- */ + +template +double MLIAPDescriptorSO3Kokkos::memory_usage() +{ + double bytes = MLIAPDescriptor::memory_usage(); + bytes += so3ptr_kokkos->memory_usage(); + + return bytes; +} + +namespace LAMMPS_NS { +template class MLIAPDescriptorSO3Kokkos; +#ifdef LMP_KOKKOS_GPU +template class MLIAPDescriptorSO3Kokkos; +#endif +} diff --git a/src/KOKKOS/mliap_descriptor_so3_kokkos.h b/src/KOKKOS/mliap_descriptor_so3_kokkos.h new file mode 100644 index 0000000000..aea4540632 --- /dev/null +++ b/src/KOKKOS/mliap_descriptor_so3_kokkos.h @@ -0,0 +1,51 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#ifndef LMP_MLIAP_DESCRIPTOR_SO3_KOKKOS_H +#define LMP_MLIAP_DESCRIPTOR_SO3_KOKKOS_H + +#include "mliap_descriptor_kokkos.h" +#include "mliap_descriptor_so3.h" +#include "mliap_so3_kokkos.h" + +namespace LAMMPS_NS { +template +class MLIAPDescriptorSO3Kokkos : + public MLIAPDescriptorSO3, + public MLIAPDescriptorKokkos { + public: + MLIAPDescriptorSO3Kokkos(LAMMPS *, char *); + ~MLIAPDescriptorSO3Kokkos() override; + + void compute_descriptors(class MLIAPData *) override; + void compute_forces(class MLIAPData *) override; + void compute_force_gradients(class MLIAPData *) override; + void compute_descriptor_gradients(class MLIAPData *) override; + void init() override; + double memory_usage() override; + + protected: + template + KOKKOS_FUNCTION static void v_tally(int vflag_either, int vflag_global, int vflag_atom, int i, + int j, int ij, double *fij, ViewType rij, + Kokkos::View virial, ViewType vatom); + class MLIAP_SO3Kokkos *so3ptr_kokkos; + + // inherited from non-Kokkos class +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/KOKKOS/mliap_model_kokkos.h b/src/KOKKOS/mliap_model_kokkos.h new file mode 100644 index 0000000000..792bad4e48 --- /dev/null +++ b/src/KOKKOS/mliap_model_kokkos.h @@ -0,0 +1,54 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#ifndef LMP_MLIAP_MODEL_KOKKOS_H +#define LMP_MLIAP_MODEL_KOKKOS_H + +#include "kokkos_type.h" +#include "memory_kokkos.h" +#include "pointers.h" + +namespace LAMMPS_NS { + +template class MLIAPModelKokkos : protected Pointers { + public: + MLIAPModelKokkos(LAMMPS *lmp, MLIAPModel *model_in) : Pointers(lmp), model(model_in) {} + virtual ~MLIAPModelKokkos() + { + memoryKK->destroy_kokkos(k_coeffelem); + model->coeffelem = nullptr; + } + + void set_k_coeffelem() + { + double **tmp = nullptr; + memoryKK->destroy_kokkos(k_coeffelem); + memoryKK->create_kokkos(k_coeffelem, tmp, model->nelements, model->nparams, + "MLIAPModelKokkos::coeffelem"); + for (int i = 0; i < model->nelements; ++i) + for (int j = 0; j < model->nparams; ++j) tmp[i][j] = model->coeffelem[i][j]; + delete model->coeffelem; + model->coeffelem = tmp; + k_coeffelem.modify(); + k_coeffelem.sync(); + } + + MLIAPModel *model; + DAT::tdual_float_2d k_coeffelem; +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/KOKKOS/mliap_model_linear_kokkos.cpp b/src/KOKKOS/mliap_model_linear_kokkos.cpp new file mode 100644 index 0000000000..14166ab3cd --- /dev/null +++ b/src/KOKKOS/mliap_model_linear_kokkos.cpp @@ -0,0 +1,99 @@ +// clang-format off +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#include "mliap_model_linear_kokkos.h" + +#include "mliap_data_kokkos.h" +#include "error.h" + +using namespace LAMMPS_NS; + +template +MLIAPModelLinearKokkos::MLIAPModelLinearKokkos(LAMMPS *lmp, char *args) : + MLIAPModelLinear(lmp,args), + MLIAPModelKokkos(lmp, this) +{ + if (args) MLIAPModelKokkos::set_k_coeffelem(); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPModelLinearKokkos::compute_gradients(class MLIAPData *data) +{ + MLIAPDataKokkos *k_data = (MLIAPDataKokkos*)(data); + + // read but never changes + auto d_coeffelem = this->k_coeffelem.template view(); + + // read + auto d_ielems = k_data->k_ielems.template view(); + auto d_descriptors = k_data->k_descriptors.template view(); + + // written + auto d_betas = k_data->k_betas.template view(); + auto d_eatoms = k_data->k_eatoms.template view(); + + const auto eflag = data->eflag; + const int ndescriptors=data->ndescriptors; + + Kokkos::parallel_reduce(data->nlistatoms, KOKKOS_LAMBDA (int ii, double &update) { + const int ielem = d_ielems(ii); + + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) + d_betas(ii,icoeff) = d_coeffelem(ielem,icoeff+1); + + // add in contributions to global and per-atom energy + // this is optional and has no effect on force calculation + if (eflag) { + // energy of atom I + double etmp = d_coeffelem(ielem,0); + + // E_i = beta.B_i + for (int icoeff = 0; icoeff < ndescriptors; icoeff++) + etmp += d_coeffelem(ielem,icoeff+1)*d_descriptors(ii, icoeff); + update += etmp; + d_eatoms(ii) = etmp; + } + }, data->energy); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPModelLinearKokkos::compute_gradgrads(class MLIAPData *data) +{ + MLIAPModelLinear::compute_gradgrads(data); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAPModelLinearKokkos::compute_force_gradients(class MLIAPData *data) +{ + MLIAPModelLinear::compute_force_gradients(data); +} + +/* ---------------------------------------------------------------------- */ + +namespace LAMMPS_NS { +template class MLIAPModelLinearKokkos; +#ifdef LMP_KOKKOS_GPU +template class MLIAPModelLinearKokkos; +#endif +} diff --git a/src/KOKKOS/mliap_model_linear_kokkos.h b/src/KOKKOS/mliap_model_linear_kokkos.h new file mode 100644 index 0000000000..829e3c4039 --- /dev/null +++ b/src/KOKKOS/mliap_model_linear_kokkos.h @@ -0,0 +1,37 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#ifndef LMP_MLIAP_MODEL_LINEAR_KOKKOS_H +#define LMP_MLIAP_MODEL_LINEAR_KOKKOS_H + +#include "mliap_model_linear.h" + +#include "mliap_model_kokkos.h" + +namespace LAMMPS_NS { + +template +class MLIAPModelLinearKokkos : public MLIAPModelLinear, public MLIAPModelKokkos { + public: + MLIAPModelLinearKokkos(LAMMPS *, char * = nullptr); + + void compute_gradients(class MLIAPData *) override; + void compute_gradgrads(class MLIAPData *) override; + void compute_force_gradients(class MLIAPData *) override; +}; +} // namespace LAMMPS_NS +#endif diff --git a/src/KOKKOS/mliap_so3_kokkos.cpp b/src/KOKKOS/mliap_so3_kokkos.cpp new file mode 100644 index 0000000000..128e0f0a0c --- /dev/null +++ b/src/KOKKOS/mliap_so3_kokkos.cpp @@ -0,0 +1,1194 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. + ------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#include "mliap_so3_kokkos.h" + +#include "error.h" +#include "math_const.h" +#include "math_special_kokkos.h" +#include "memory.h" +#include "memory_kokkos.h" +#include "mliap_so3_math.h" + +#include + +using namespace SO3Math; +using namespace LAMMPS_NS; +using namespace MathConst; +using namespace MathSpecialKokkos; + +#define SMALL 1.0e-8 + +/* ---------------------------------------------------------------------- */ + +template +MLIAP_SO3Kokkos::MLIAP_SO3Kokkos(LAMMPS *lmp, double vrcut, int vlmax, int vnmax, double valpha) : Pointers(lmp) +{ + m_rcut = vrcut; + m_alpha = valpha; + m_lmax = vlmax; + m_nmax = vnmax; + compute_ncoeff(); + + m_Nmax = (m_nmax + m_lmax + 1) * 10; + m_numYlms = (m_lmax + 1) * (m_lmax + 1); + + m_init_arrays = 0; + m_dfac_l1 = m_dfac_l2 = 0; + m_pfac_l1 = m_pfac_l2 = 0; + m_idxu_count = m_idxy_count = 0; + alloc_init = alloc_arrays = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +template +MLIAP_SO3Kokkos::~MLIAP_SO3Kokkos() +{ + memoryKK->destroy_kokkos(m_ellpl1); + memoryKK->destroy_kokkos(m_ellm1); + memoryKK->destroy_kokkos(m_pfac); + memoryKK->destroy_kokkos(m_Ylms); + memoryKK->destroy_kokkos(m_dfac0); + memoryKK->destroy_kokkos(m_dfac1); + memoryKK->destroy_kokkos(m_dfac2); + memoryKK->destroy_kokkos(m_dfac3); + memoryKK->destroy_kokkos(m_dfac4); + memoryKK->destroy_kokkos(m_dfac5); + memoryKK->destroy_kokkos(m_w); + memoryKK->destroy_kokkos(m_g_array); + + memoryKK->destroy_kokkos(m_rootpq); + memoryKK->destroy_kokkos(m_idxu_block); + memoryKK->destroy_kokkos(m_idxylm); + + memoryKK->destroy_kokkos(m_rip_array); + memoryKK->destroy_kokkos(m_rip_darray); + + memoryKK->destroy_kokkos(m_sbes_array); + memoryKK->destroy_kokkos(m_sbes_darray); + + memoryKK->destroy_kokkos(m_plist_r); + + memoryKK->destroy_kokkos(m_ulist_r); + memoryKK->destroy_kokkos(m_ulist_i); + + memoryKK->destroy_kokkos(m_dYlm_r); + memoryKK->destroy_kokkos(m_dYlm_i); + + memoryKK->destroy_kokkos(k_dplist_r); + + memoryKK->destroy_kokkos(m_dclist); + + memoryKK->destroy_kokkos(m_clisttot_r); + memoryKK->destroy_kokkos(m_clisttot_i); + + t_numneighs = int_1d(); + t_jelems = int_1d(); + t_wjelem = float_1d(); + t_rij = float_2d(); + t_ij = int_1d(); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::compute_ncoeff() +{ + ncoeff = m_nmax * (m_nmax + 1) * (m_lmax + 1) / 2; +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::init() +{ + int totali; + + totali = m_lmax + 1; + memoryKK->destroy_kokkos(m_ellpl1); + memoryKK->create_kokkos(m_ellpl1, totali, "MLIAP_SO3Kokkos:m_ellpl1"); + memoryKK->destroy_kokkos(m_ellm1); + memoryKK->create_kokkos(m_ellm1, totali, "MLIAP_SO3Kokkos:m_ellm1"); + alloc_init = 2.0 * totali * sizeof(double); + using range=Kokkos::RangePolicy; + auto ellpl1 = m_ellpl1, ellm1 = m_ellm1; + + Kokkos::parallel_for(range(0,m_lmax), KOKKOS_LAMBDA (int ll) { + int l=ll+1; + ellpl1[l] = get_sum(0, l + 2, 1, 2); + ellm1[l] = get_sum(0, l, 1, 2); + }); + + double pfac1 = 1.0 / 4.0 / MY_PI; + m_pfac_l1 = m_lmax + 2; + m_pfac_l2 = (m_lmax + 2) * (m_lmax + 2) + 1; + totali = m_pfac_l1 * m_pfac_l2; + memoryKK->destroy_kokkos(m_pfac); + memoryKK->create_kokkos(m_pfac, totali, "MLIAP_SO3Kokkos:m_pfac"); + memoryKK->destroy_kokkos(m_Ylms); + memoryKK->create_kokkos(m_Ylms, totali, "MLIAP_SO3Kokkos:m_Ylms"); + alloc_init += 2 * totali * sizeof(double); + + auto pfac = m_pfac, Ylms = m_Ylms; + auto pfac_l2 = m_pfac_l2, lmax = m_lmax; + // Serial but just to make sure run with device memory + Kokkos::parallel_for(range(0,1), KOKKOS_LAMBDA (int ) { + int i=0; + for (int l = 0; l < lmax + 2; l++) + for (int m = -l; m < l + 1; m++) { + pfac[l * pfac_l2 + m] = sqrt((2.0 * l + 1.0) * pfac1) * powsign(m); + Ylms[i] = pfac[l * pfac_l2 + m]; + i += 1; + } + }); + + m_dfac_l1 = m_lmax + 1; + m_dfac_l2 = m_numYlms + 1; + totali = m_dfac_l1 * m_dfac_l2; + memoryKK->destroy_kokkos(m_dfac0); + memoryKK->create_kokkos(m_dfac0, totali, "MLIAP_SO3Kokkos:m_dfac0"); + memoryKK->destroy_kokkos(m_dfac1); + memoryKK->create_kokkos(m_dfac1, totali, "MLIAP_SO3Kokkos:m_dfac1"); + memoryKK->destroy_kokkos(m_dfac2); + memoryKK->create_kokkos(m_dfac2, totali, "MLIAP_SO3Kokkos:m_dfac2"); + memoryKK->destroy_kokkos(m_dfac3); + memoryKK->create_kokkos(m_dfac3, totali, "MLIAP_SO3Kokkos:m_dfac3"); + memoryKK->destroy_kokkos(m_dfac4); + memoryKK->create_kokkos(m_dfac4, totali, "MLIAP_SO3Kokkos:m_dfac4"); + memoryKK->destroy_kokkos(m_dfac5); + memoryKK->create_kokkos(m_dfac5, totali, "MLIAP_SO3Kokkos:m_dfac5"); + alloc_init += 6.0 * totali * sizeof(double); + + auto dfac0 = m_dfac0,dfac1 = m_dfac1,dfac2 = m_dfac2,dfac3 = m_dfac3,dfac4 = m_dfac4,dfac5 = m_dfac5; + auto dfac_l2 = m_dfac_l2; + + Kokkos::parallel_for(range(0,m_lmax), KOKKOS_LAMBDA (int ll) { + int l = ll+1; + for (int m = -l; m < l + 1; m++) { + dfac0[l * dfac_l2 + m] = + -sqrt(((l + 1.0) * (l + 1.0) - m * m) / (2.0 * l + 1.0) / (2.0 * l + 3.0)) * l; + dfac1[l * dfac_l2 + m] = + sqrt((l * l - m * m) / (2.0 * l - 1.0) / (2.0 * l + 1.0)) * (l + 1.0); + dfac2[l * dfac_l2 + m] = + -sqrt((l + m + 1.0) * (l + m + 2.0) / 2.0 / (2.0 * l + 1.0) / (2.0 * l + 3.0)) * l; + dfac3[l * dfac_l2 + m] = + sqrt((l - m - 1.0) * (l - m) / 2.0 / (2.0 * l - 1.0) / (2.0 * l + 1.0)) * (l + 1.0); + dfac4[l * dfac_l2 + m] = + -sqrt((l - m + 1.0) * (l - m + 2.0) / 2.0 / (2.0 * l + 1.0) / (2.0 * l + 3.0)) * l; + dfac5[l * dfac_l2 + m] = + sqrt((l + m - 1.0) * (l + m) / 2.0 / (2.0 * l - 1.0) / (2.0 * l + 1.0)) * (l + 1.0); + } + }); + + totali = m_nmax * m_nmax; + memoryKK->destroy_kokkos(m_w); + memoryKK->create_kokkos(m_w, totali, "MLIAP_SO3Kokkos:w"); + alloc_init += totali * sizeof(double); + + totali = m_nmax * m_Nmax; + memoryKK->destroy_kokkos(m_g_array); + memoryKK->create_kokkos(m_g_array, totali, "MLIAP_SO3Kokkos:g_array"); + alloc_init += totali * sizeof(double); + + { + Kokkos::View w("w", m_nmax * m_nmax), g_array("g_array", m_nmax * m_Nmax); + compute_W(m_nmax, w.data()); + init_garray(m_nmax, m_lmax, m_rcut, m_alpha, w.data(), m_nmax, g_array.data(), m_Nmax); + Kokkos::deep_copy(m_w,w); + Kokkos::deep_copy(m_g_array,g_array); + } + + int twolmax; + twolmax = 2 * (m_lmax + 1); + int m_ldim = twolmax + 1; + totali = m_ldim * m_ldim; + memoryKK->destroy_kokkos(m_rootpq); + memoryKK->create_kokkos(m_rootpq, totali, "MLIAP_SO3Kokkos:rootpq"); + alloc_init += totali * sizeof(double); + + auto rootpq=m_rootpq; + auto ldim=m_ldim; + Kokkos::parallel_for(range(1,m_ldim), KOKKOS_LAMBDA (int p) { + for (int q = 1; q < ldim; q++) + rootpq[p * ldim + q] = sqrt(static_cast(p) / q); + }); + + memoryKK->destroy_kokkos(m_idxu_block); + memoryKK->create_kokkos(m_idxu_block, m_ldim, "MLIAP_SO3Kokkos:idxu_bloc"); + alloc_init += totali * sizeof(double); + + totali = square(m_lmax + 2); + memoryKK->destroy_kokkos(m_idxylm); + memoryKK->create_kokkos(m_idxylm, totali, "MLIAP_SO3Kokkos:idxylm"); + alloc_init += totali * sizeof(double); + + Kokkos::View idxu_block("idxu_block",m_ldim); + Kokkos::View idxylm("idxylm",totali); + m_idxu_count = m_idxy_count = 0; + for (int l = 0; l < m_ldim; l++) { + idxu_block[l] = m_idxu_count; + for (int mb = 0; mb < l + 1; mb++) + for (int ma = 0; ma < l + 1; ma++) { + if (l % 2 == 0 && ma == l / 2) { + idxylm[m_idxy_count] = m_idxu_count; + m_idxy_count += 1; + } + m_idxu_count += 1; + } + } + Kokkos::deep_copy(m_idxu_block,idxu_block); + Kokkos::deep_copy(m_idxylm,idxylm); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::init_arrays(int nlocal, int ncoefs) +{ + + { + // sum up the temp memory space per particle + int sum_of_temps=0; + //ulist_r and i + sum_of_temps += m_idxu_count*2; + //m_dYlm_r + sum_of_temps += m_numYlms * 3 * 2; + //m_clisttot_r + sum_of_temps += m_nmax * m_numYlms * 2; + // k_dclist + sum_of_temps += m_nmax * m_numYlms * 3 * 2; + m_chunk_size = m_temp_memory_size/(sum_of_temps * sizeof(double)); + } + + int totali = nlocal * ncoefs; + if ( nlocal > (int)m_plist_r.extent(0)) { + memoryKK->destroy_kokkos(m_plist_r); + memoryKK->create_kokkos(m_plist_r, nlocal, ncoefs, "MLIAP_SO3Kokkos:m_plist_r"); + alloc_arrays = totali * sizeof(double); + } + + int num_of_temp = std::min(nlocal, m_chunk_size); + if ((int)m_ulist_r.extent(0) < num_of_temp ) { + totali = m_idxu_count; + memoryKK->destroy_kokkos(m_ulist_r); + memoryKK->create_kokkos(m_ulist_r, num_of_temp, totali, "MLIAP_SO3Kokkos:m_ulist_r"); + memoryKK->destroy_kokkos(m_ulist_i); + memoryKK->create_kokkos(m_ulist_i, num_of_temp, totali, "MLIAP_SO3Kokkos:m_ulist_i"); + alloc_arrays += 2.0 * totali * num_of_temp * sizeof(double); + + totali = m_numYlms * 3; + memoryKK->destroy_kokkos(m_dYlm_r); + memoryKK->create_kokkos(m_dYlm_r, num_of_temp, m_numYlms, 3, "MLIAP_SO3Kokkos:m_dYlm_r"); + memoryKK->destroy_kokkos(m_dYlm_i); + memoryKK->create_kokkos(m_dYlm_i, num_of_temp, m_numYlms, 3, "MLIAP_SO3Kokkos:m_dYlm_i"); + alloc_arrays += 2.0 * m_numYlms * 3 * num_of_temp * sizeof(double); + + memoryKK->destroy_kokkos(m_dclist); + memoryKK->create_kokkos(m_dclist, num_of_temp, m_nmax, m_numYlms, 3, "MLIAP_SO3Kokkos:k_dclist_r"); + alloc_arrays += m_nmax * m_numYlms * 3 * num_of_temp* sizeof(double); + + memoryKK->destroy_kokkos(m_clisttot_r); + memoryKK->create_kokkos(m_clisttot_r, num_of_temp, m_nmax, m_numYlms, "MLIAP_SO3Kokkos:m_clisttot_r"); + memoryKK->destroy_kokkos(m_clisttot_i); + memoryKK->create_kokkos(m_clisttot_i, num_of_temp, m_nmax, m_numYlms, "MLIAP_SO3Kokkos:m_clisttot_i"); + alloc_arrays += 2.0 * m_nmax * m_numYlms * num_of_temp * sizeof(double); + m_init_arrays = 1; + } +} + +/* ---------------------------------------------------------------------- */ + +template +double MLIAP_SO3Kokkos::memory_usage() +{ + return alloc_init + alloc_arrays; +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::compute_W(int nmax, double *arr) +{ + int alpha, beta, temp1, temp2; + + for (alpha = 1; alpha < nmax + 1; alpha++) { + temp1 = (2 * alpha + 5) * (2 * alpha + 6) * (2 * alpha + 7); + for (beta = 1; beta < alpha + 1; beta++) { + temp2 = (2 * beta + 5) * (2 * beta + 6) * (2 * beta + 7); + arr[(alpha - 1) * nmax + beta - 1] = + sqrt(temp1 * temp2) / (5 + alpha + beta) / (6 + alpha + beta) / (7 + alpha + beta); + arr[(beta - 1) * nmax + alpha - 1] = arr[(alpha - 1) * nmax + beta - 1]; + } + } + + int i, j, k, n = nmax; + auto outeig = new double[n]; + auto outeigvec = new double[n * n]; + auto arrinv = new double[n * n]; + + auto sqrtD = new double[n * n]; + auto tempM = new double[n * n]; + + auto temparr = new double *[n]; + auto tempvl = new double *[n]; + auto tempout = new double[n]; + + int info; + + info = invert_matrix(n, arr, arrinv); + if (info != 0) error->all(FLERR, "Invert matrix Error in W calculation!"); + + for (i = 0; i < n; i++) { + temparr[i] = new double[n]; + tempvl[i] = new double[n]; + for (j = 0; j < n; j++) temparr[i][j] = arrinv[i * n + j]; + } + + jacobin(n, temparr, tempout, tempvl); + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) outeigvec[i * n + j] = tempvl[i][j]; + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) + if (i == j) + sqrtD[i * n + j] = sqrt(tempout[i]); + else + sqrtD[i * n + j] = 0.0; + + double dtemp; + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) { + dtemp = 0; + for (k = 0; k < n; k++) dtemp += outeigvec[i * n + k] * sqrtD[k * n + j]; + + tempM[i * n + j] = dtemp; + } + + info = invert_matrix(n, outeigvec, arrinv); + if (info != 0) error->all(FLERR, "Invert matrix Error in W calculation!"); + + for (i = 0; i < n; i++) + for (j = 0; j < n; j++) { + dtemp = 0; + for (k = 0; k < n; k++) dtemp += tempM[i * n + k] * arrinv[k * n + j]; + + arr[i * n + j] = dtemp; + } + + delete[] outeig; + delete[] outeigvec; + delete[] arrinv; + + delete[] sqrtD; + delete[] tempM; + + for (i = 0; i < n; i++) { + delete[] temparr[i]; + delete[] tempvl[i]; + } + delete[] temparr; + delete[] tempvl; + delete[] tempout; +} + +/* ---------------------------------------------------------------------- */ +template +template +KOKKOS_INLINE_FUNCTION +void MLIAP_SO3Kokkos::compute_pi(int nmax, int lmax, ViewType clisttot_r, ViewType clisttot_i, int /*lcl2*/, + float_2d plist_r, int indpl) const +{ + int n1, n2, j, l, m, i = 0; + double norm; + for (n1 = 0; n1 < nmax; n1++) + for (n2 = 0; n2 < n1 + 1; n2++) { + j = 0; + for (l = 0; l < lmax + 1; l++) { + norm = 2.0 * sqrt(2.0) * MY_PI / sqrt(2.0 * l + 1.0); + + for (m = -l; m < l + 1; m++) { + + plist_r(indpl, i) += (clisttot_r(n1, j) * clisttot_r(n2, j) + + clisttot_i(n1, j) * clisttot_i(n2, j)) * + norm; + j += 1; + } + i += 1; + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +double MLIAP_SO3Kokkos::phi(double r, int alpha, double rcut) +{ + return powint((rcut - r), (alpha + 2)) / + sqrt(2 * powint(rcut, (2 * alpha + 7)) / (2 * alpha + 5) / (2 * alpha + 6) / (2 * alpha + 7)); +} + +/* ---------------------------------------------------------------------- */ + +template +double MLIAP_SO3Kokkos::compute_g(double r, int n, int nmax, double rcut, double *w, int lw1) +{ + double Sum; + Sum = 0.0; + int alpha; + + for (alpha = 1; alpha < nmax + 1; alpha++) + Sum += w[(n - 1) * lw1 + alpha - 1] * phi(r, alpha, rcut); + + return Sum; +} + +/* ---------------------------------------------------------------------- */ +template +KOKKOS_INLINE_FUNCTION +double MLIAP_SO3Kokkos::Cosine(double Rij, double Rc) const +{ + + return 0.5 * (cos(MY_PI * Rij / Rc) + 1.0); +} + +/* ---------------------------------------------------------------------- */ +template +KOKKOS_INLINE_FUNCTION +double MLIAP_SO3Kokkos::CosinePrime(double Rij, double Rc) const +{ + + return -0.5 * MY_PI / Rc * sin(MY_PI * Rij / Rc); +} + +/* ---------------------------------------------------------------------- */ +template +KOKKOS_INLINE_FUNCTION +double MLIAP_SO3Kokkos::compute_sfac(double r, double rcut) const +{ + if (r > rcut) + return 0.0; + else + return Cosine(r, rcut); +} + +/* ---------------------------------------------------------------------- */ +template +KOKKOS_INLINE_FUNCTION +double MLIAP_SO3Kokkos::compute_dsfac(double r, double rcut) const +{ + if (r > rcut) + return 0.0; + else + return CosinePrime(r, rcut); +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +int MLIAP_SO3Kokkos::get_sum(int istart, int iend, int id, int imult) +{ + int ires = 0; + int i; + + for (i = istart; i < iend; i = i + id) ires += i * imult; + + return ires; +} + +/* ---------------------------------------------------------------------- */ + +template +template +KOKKOS_INLINE_FUNCTION +void MLIAP_SO3Kokkos::compute_uarray_recursive(double x, double y, double z, double r, int twol, + UlistView ulist_r, UlistView ulist_i, int_1d idxu_block, + float_1d rootpqarray) const +{ + int l, llu, llup, mb, ma, mbpar, mapar; + double rootpq; + int ldim = twol + 1; + + double theta, phi, atheta, btheta; + + double aphi_r, aphi_i, a_r, a_i, b_r, b_i; + + theta = acos(z / r); + phi = atan2(y, x); + + atheta = cos(theta / 2); + btheta = sin(theta / 2); + + aphi_r = cos(phi / 2); + aphi_i = sin(phi / 2); + + a_r = atheta * aphi_r; + a_i = atheta * aphi_i; + b_r = btheta * aphi_r; + b_i = btheta * aphi_i; + + ulist_r[0] = 1.0; + ulist_i[0] = 0.0; + + for (l = 1; l < ldim; l++) { + + llu = idxu_block[l]; + llup = idxu_block[l - 1]; + mb = 0; + + while (2 * mb <= l) { + + ulist_r[llu] = 0.0; + ulist_i[llu] = 0.0; + for (ma = 0; ma < l; ma++) { + + rootpq = rootpqarray[(l - ma) * ldim + l - mb]; + + ulist_r[llu] += rootpq * (a_r * ulist_r[llup] + a_i * ulist_i[llup]); + ulist_i[llu] += rootpq * (a_r * ulist_i[llup] - a_i * ulist_r[llup]); + + rootpq = rootpqarray[(ma + 1) * ldim + l - mb]; + + ulist_r[llu + 1] += -rootpq * (b_r * ulist_r[llup] + b_i * ulist_i[llup]); + ulist_i[llu + 1] += -rootpq * (b_r * ulist_i[llup] - b_i * ulist_r[llup]); + + llu += 1; + llup += 1; + } + + llu += 1; + mb += 1; + } + + llu = idxu_block[l]; + llup = llu + (l + 1) * (l + 1) - 1; + mbpar = 1; + mb = 0; + + while (2 * mb <= l) { + mapar = mbpar; + for (ma = 0; ma < l + 1; ma++) { + if (mapar == 1) { + + ulist_r[llup] = ulist_r[llu]; + ulist_i[llup] = -ulist_i[llu]; + + } else { + + ulist_r[llup] = -ulist_r[llu]; + ulist_i[llup] = ulist_i[llu]; + } + mapar = -mapar; + llu += 1; + llup -= 1; + } + mbpar = -mbpar; + mb += 1; + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::init_garray(int nmax, int lmax, double rcut, double alpha, double *w, int lw1, + double *g_array, int lg2) +{ + int i, n, Nmax = (nmax + lmax + 1) * 10; + double x, xi; + + for (i = 1; i < Nmax + 1; i++) { + // roots of Chebyshev polynomial of degree N + x = cos((2 * i - 1) * MY_PI / 2 / Nmax); + // transform the interval [-1,1] to [0, rcut] + xi = rcut / 2 * (x + 1); + for (n = 1; n < nmax + 1; n++) + // r**2*g(n)(r)*e^(-alpha*r**2) + g_array[(n - 1) * lg2 + i - 1] = rcut / 2 * MY_PI / Nmax * sqrt(1 - x * x) * xi * xi * + compute_g(xi, n, nmax, rcut, w, lw1) * exp(-alpha * xi * xi); + } +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void MLIAP_SO3Kokkos::operator() (const MLIAPSO3GetSBESArrayTag&, int ii) const{ + int ipair = t_ij(ii); + for (int neighbor = 0; neighbor < t_numneighs[ii]; neighbor++) { + double x = t_rij(ipair, 0); + double y = t_rij(ipair, 1); + double z = t_rij(ipair, 2); + ipair++; + + double ri = sqrt(x * x + y * y + z * z); + + if (ri < SMALL) continue; + + const double pfac1 = t_alpha * t_rcut; + const double pfac4 = t_rcut / 2.0; + const double pfac3 = MY_PI / 2.0 / m_Nmax; + double pfac2 = pfac1 * ri; + const int findex = m_Nmax * (m_lmax + 1); + const int gindex = (ipair - 1) * findex; + const int mindex = m_lmax + 1; + + for (int i = 1; i < m_Nmax + 1; i++) { + const bigint i1mindex = (bigint) (i - 1) * mindex; + + x = cos((2 * i - 1) * pfac3); + double xi = pfac4 * (x + 1); + double rb = pfac2 * (x + 1); + + double sa = sinh(rb) / rb; + double sb = (cosh(rb) - sa) / rb; + + m_sbes_array[gindex + i1mindex + 0] = sa; + m_sbes_array[gindex + i1mindex + 1] = sb; + + int j; + for (j = 2; j < t_lmax + 1; j++) + m_sbes_array[gindex + i1mindex + j] = m_sbes_array[gindex + i1mindex + j - 2] - + (2 * j - 1) / rb * m_sbes_array[gindex + i1mindex + j - 1]; + + double exts = m_sbes_array[gindex + i1mindex + j - 2] - + (2 * j - 1) / rb * m_sbes_array[gindex + i1mindex + j - 1]; + + m_sbes_darray[gindex + i1mindex + 0] = sb; + + for (j = 1; j < t_lmax; j++) + m_sbes_darray[gindex + i1mindex + j] = xi * + (j * m_sbes_array[gindex + i1mindex + j - 1] + + (j + 1) * m_sbes_array[gindex + i1mindex + j + 1]) / + (2 * j + 1); + + m_sbes_darray[gindex + i1mindex + j] = + xi * (j * m_sbes_array[gindex + i1mindex + j - 1] + (j + 1) * exts) / (2 * j + 1); + m_sbes_darray[gindex + i1mindex + 0] = xi * sb; + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void MLIAP_SO3Kokkos::operator() (const MLIAPSO3GetRipArrayTag&, int ii) const{ + int ipair = t_ij(ii); + for (int neighbor = 0; neighbor < t_numneighs[ii]; neighbor++) { + + double x = t_rij(ipair, 0); + double y = t_rij(ipair, 1); + double z = t_rij(ipair, 2); + ipair++; + + double ri = sqrt(x * x + y * y + z * z); + if (ri < SMALL) continue; + + double expfac = 4 * MY_PI * exp(-t_alpha * ri * ri); + for (int n = 1; n < t_nmax + 1; n++) + for (int l = 0; l < t_lmax + 1; l++) { + double integral = 0.0, integrald = 0.0; + for (int i = 0; i < m_Nmax; i++) { + integral += m_g_array[(n - 1) * m_Nmax + i] * + m_sbes_array[(ipair - 1) * m_Nmax * (m_lmax + 1) + (bigint) i * (m_lmax + 1) + l]; + integrald += m_g_array[(n - 1) * m_Nmax + i] * + m_sbes_darray[(ipair - 1) * m_Nmax * (m_lmax + 1) + (bigint) i * (m_lmax + 1) + l]; + } + + m_rip_array[(ipair - 1) * m_nmax * (m_lmax + 1) + (bigint) (n - 1) * (m_lmax + 1) + l] = + integral * expfac; + m_rip_darray[(ipair - 1) * m_nmax * (m_lmax + 1) + (bigint) (n - 1) * (m_lmax + 1) + l] = + integrald * expfac; + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::spectrum(int nlocal, DAT::tdual_int_1d numneighs, DAT::tdual_int_1d jelems, DAT::tdual_float_1d wjelem, + DAT::tdual_float_2d rij, DAT::tdual_int_1d k_ij, + int nmax, int lmax, double rcut, double alpha, int totaln, int ncoefs) +{ + init_arrays(nlocal, ncoefs); + + bigint totali; + + totali = totaln * m_Nmax * (m_lmax + 1); + if ( totali > (int)m_sbes_array.extent(0)) { + memoryKK->realloc_kokkos(m_sbes_array, "MLIAP_SO3Kokkos:m_sbes_array", totali); + memoryKK->realloc_kokkos(m_sbes_darray, "MLIAP_SO3Kokkos:m_sbes_darray", totali); + alloc_arrays += 2.0 * totali * sizeof(double); + } + + totali = totaln * m_nmax * (m_lmax + 1); + if ( totali > (int)m_rip_array.extent(0)) { + memoryKK->realloc_kokkos(m_rip_array, "MLIAP_SO3Kokkos:m_rip_array", totali); + memoryKK->realloc_kokkos(m_rip_darray, "MLIAP_SO3Kokkos:m_rip_darray", totali); + alloc_arrays += 2.0 * totali * sizeof(double); + } + + totali = totaln * ncoefs * 3; + if ( totali > (int)k_dplist_r.extent(0)) { + memoryKK->realloc_kokkos(k_dplist_r, "MLIAP_SO3Kokkos:m_dplist_r", (int)totaln, ncoefs, 3); + alloc_arrays += 2.0 * totali * sizeof(double); + } + + t_numneighs = numneighs.template view(); + t_jelems = jelems.template view(); + t_wjelem = wjelem.template view(); + t_rij = rij.template view(); + t_ij = k_ij.template view(); + t_nmax = nmax; + t_lmax = lmax; + t_rcut = rcut; + t_alpha = alpha; + + { + Kokkos::RangePolicy range(0,nlocal); + Kokkos::parallel_for(range, *this); + } + + { + Kokkos::RangePolicy range(0,nlocal); + Kokkos::parallel_for(range, *this); + } + + Kokkos::deep_copy(m_plist_r,0.); + for (int start=0; start < nlocal; start += m_chunk_size) { + int end=std::min(nlocal, start+m_chunk_size); + Kokkos::deep_copy(m_clisttot_r, 0); + Kokkos::deep_copy(m_clisttot_i, 0); + { + Kokkos::RangePolicy range(start,end); + Kokkos::parallel_for(range, *this); + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void MLIAP_SO3Kokkos::operator() (const MLIAP_SO3Kokkos::MLIAPSO3SpectrumTag&, int ii) const { + int ii_chunk = ii%m_chunk_size; + auto ulist_r = Kokkos::subview(m_ulist_r,ii_chunk,Kokkos::ALL); + auto ulist_i = Kokkos::subview(m_ulist_i,ii_chunk,Kokkos::ALL); + auto clisttot_r=Kokkos::subview(m_clisttot_r,ii_chunk,Kokkos::ALL, Kokkos::ALL); + auto clisttot_i=Kokkos::subview(m_clisttot_i,ii_chunk,Kokkos::ALL, Kokkos::ALL); + + const int twolmax = 2 * (t_lmax + 1); + const int findex = m_nmax * (m_lmax + 1); + int ipair = t_ij(ii); + for (int neighbor = 0; neighbor < t_numneighs[ii]; neighbor++) { + const int jelem = t_jelems[ipair]; + int weight = t_wjelem[jelem]; + + double x = t_rij(ipair, 0); + double y = t_rij(ipair, 1); + double z = t_rij(ipair, 2); + ipair++; + + double r = sqrt(x * x + y * y + z * z); + + if (r < SMALL) continue; + for (int ti = 0; ti < m_idxu_count; ti++) { + ulist_r[ti] = 0.0; + ulist_i[ti] = 0.0; + } + + compute_uarray_recursive(x, y, z, r, twolmax, ulist_r, ulist_i, m_idxu_block, m_rootpq); + + double sfac_weight = compute_sfac(r, t_rcut)*double(weight); + + int gindex = (ipair - 1) * findex; + for (int n = 1; n < t_nmax + 1; n++) { + int i = 0; + for (int l = 0; l < t_lmax + 1; l++) { + double r_int = m_rip_array[gindex + (bigint) (n - 1) * (m_lmax + 1) + l]; + + for (int m = -l; m < l + 1; m++) { + + double Ylm_r = (ulist_r[m_idxylm[i]]) * m_pfac[l * m_pfac_l2 + m]; + clisttot_r((n - 1), i) += r_int * Ylm_r * sfac_weight; + double Ylm_i = (ulist_i[m_idxylm[i]]) * m_pfac[l * m_pfac_l2 + m]; + clisttot_i((n - 1), i) += r_int * Ylm_i * sfac_weight; + i += 1; + } + } + } + } + compute_pi(t_nmax, t_lmax, clisttot_r, clisttot_i, m_numYlms, m_plist_r, ii); +} + +/* ---------------------------------------------------------------------- */ + +template +void MLIAP_SO3Kokkos::spectrum_dxdr(int nlocal, DAT::tdual_int_1d numneighs, DAT::tdual_int_1d jelems, DAT::tdual_float_1d wjelem, + DAT::tdual_float_2d rij, DAT::tdual_int_1d k_ij, + int nmax, int lmax, double rcut, double alpha, bigint totaln, + int ncoefs) +{ + bigint totali; + + if ( nlocal > (int)m_clisttot_r.extent(0)){ + memoryKK->destroy_kokkos(m_clisttot_r); + memoryKK->create_kokkos(m_clisttot_r, nlocal, m_nmax, m_numYlms, "MLIAP_SO3Kokkos:m_clisttot_r"); + memoryKK->destroy_kokkos(m_clisttot_i); + memoryKK->create_kokkos(m_clisttot_i, nlocal, m_nmax, m_numYlms, "MLIAP_SO3Kokkos:m_clisttot_i"); + int num_of_temp = std::min(nlocal, m_chunk_size); + int delta=num_of_temp-m_ulist_r.extent(0); + if (delta > 0 ){ + memoryKK->destroy_kokkos(m_ulist_r); + memoryKK->create_kokkos(m_ulist_r, num_of_temp, m_idxu_count, "MLIAP_SO3Kokkos:m_ulist_r"); + memoryKK->destroy_kokkos(m_ulist_i); + memoryKK->create_kokkos(m_ulist_i, num_of_temp, m_idxu_count, "MLIAP_SO3Kokkos:m_ulist_i"); + alloc_arrays += 2.0 * m_idxu_count * delta * sizeof(double); + memoryKK->destroy_kokkos(m_dYlm_r); + memoryKK->create_kokkos(m_dYlm_r, num_of_temp, m_numYlms, 3, "MLIAP_SO3Kokkos:m_dYlm_r"); + memoryKK->destroy_kokkos(m_dYlm_i); + memoryKK->create_kokkos(m_dYlm_i, num_of_temp, m_numYlms, 3, "MLIAP_SO3Kokkos:m_dYlm_i"); + alloc_arrays += 2.0 * m_numYlms * 3 * delta * sizeof(double); + } + } + + totali = totaln * m_Nmax * (m_lmax + 1); + if ( totali > (int)m_sbes_array.extent(0)) { + memoryKK->destroy_kokkos(m_sbes_array); + memoryKK->create_kokkos(m_sbes_array, totali, "MLIAP_SO3Kokkos:m_sbes_array"); + memoryKK->destroy_kokkos(m_sbes_darray); + memoryKK->create_kokkos(m_sbes_darray, totali, "MLIAP_SO3Kokkos:m_sbes_darray"); + + totali = totaln * m_nmax * (m_lmax + 1); + memoryKK->destroy_kokkos(m_rip_array); + memoryKK->create_kokkos(m_rip_array, totali, "MLIAP_SO3Kokkos:m_rip_array"); + memoryKK->destroy_kokkos(m_rip_darray); + memoryKK->create_kokkos(m_rip_darray, totali, "MLIAP_SO3Kokkos:m_rip_darray"); + + memoryKK->destroy_kokkos(k_dplist_r); + memoryKK->create_kokkos(k_dplist_r, (int)totaln, ncoefs, 3, "MLIAP_SO3Kokkos:m_dplist_r"); + } + + t_numneighs=numneighs.template view(); + t_jelems=jelems.template view(); + t_wjelem=wjelem.template view(); + t_rij=rij.template view(); + t_ij = k_ij.template view(); + t_nmax=nmax; + t_lmax=lmax; + t_rcut=rcut; + t_alpha=alpha; + + { + Kokkos::RangePolicy range(0,nlocal); + Kokkos::parallel_for(range, *this); + } + + { + Kokkos::RangePolicy range(0,nlocal); + Kokkos::parallel_for(range, *this); + } + + Kokkos::deep_copy(k_dplist_r, 0.0); + for (int start=0; start < nlocal; start += m_chunk_size) { + int end=std::min(nlocal, start+m_chunk_size); + Kokkos::deep_copy(m_clisttot_r, 0); + Kokkos::deep_copy(m_clisttot_i, 0); + { + Kokkos::RangePolicy range(start,end); + Kokkos::parallel_for(range, *this); + } + } +} + +/* ---------------------------------------------------------------------- */ + +template +KOKKOS_INLINE_FUNCTION +void MLIAP_SO3Kokkos::operator() (const MLIAP_SO3Kokkos::MLIAPSO3SpectrumDXDRTag&, int ii) const { + //TO-DO Need to move m_ulist_r, m_ulist_i into local shared memory + int ii_chunk = ii%m_chunk_size; + auto ulist_r = Kokkos::subview(m_ulist_r,ii_chunk,Kokkos::ALL); + auto ulist_i = Kokkos::subview(m_ulist_i,ii_chunk,Kokkos::ALL); + auto clisttot_r=Kokkos::subview(m_clisttot_r,ii_chunk,Kokkos::ALL, Kokkos::ALL); + auto clisttot_i=Kokkos::subview(m_clisttot_i,ii_chunk,Kokkos::ALL, Kokkos::ALL); + auto dYlm_r = Kokkos::subview(m_dYlm_r,ii_chunk,Kokkos::ALL,Kokkos::ALL); + auto dYlm_i = Kokkos::subview(m_dYlm_i,ii_chunk,Kokkos::ALL,Kokkos::ALL); + auto dclist = Kokkos::subview(m_dclist,ii_chunk,Kokkos::ALL,Kokkos::ALL,Kokkos::ALL); + int twolmax = 2 * (t_lmax + 1); + int findex = m_nmax * (m_lmax + 1); + int ipair = t_ij(ii); + for (int neighbor = 0; neighbor < t_numneighs[ii]; neighbor++) { + + const int jelem = t_jelems[ipair]; + int weight = t_wjelem[jelem]; + + double x = t_rij(ipair, 0); + double y = t_rij(ipair, 1); + double z = t_rij(ipair, 2); + ipair++; + + double r = sqrt(x * x + y * y + z * z); + + if (r < SMALL) continue; + + for (bigint ti = 0; ti < m_idxu_count; ti++) { + ulist_r[ti] = 0.0; + ulist_i[ti] = 0.0; + } + + compute_uarray_recursive(x, y, z, r, twolmax, ulist_r, ulist_i, m_idxu_block, m_rootpq); + + double sfac_weight = compute_sfac(r, t_rcut)*double(weight); + bigint gindex = (ipair - 1) * findex; + for (int n = 0; n < t_nmax; n++) { + int i = 0; + for (int l = 0; l < t_lmax + 1; l++) { + double r_int = m_rip_array[gindex + n * (m_lmax + 1) + l]; + + for (int m = -l; m < l + 1; m++) { + double Ylm_r = (ulist_r[m_idxylm[i]]) * m_pfac[l * m_pfac_l2 + m]; + clisttot_r(n, i) += r_int * Ylm_r * sfac_weight; + double Ylm_i = (ulist_i[m_idxylm[i]]) * m_pfac[l * m_pfac_l2 + m]; + clisttot_i(n, i) += r_int * Ylm_i * sfac_weight; + i += 1; + } + } + } + + } + + ipair = t_ij(ii); + for (int neighbor = 0; neighbor < t_numneighs[ii]; neighbor++) { + const int jelem = t_jelems[ipair]; + int weight = t_wjelem[jelem]; + + double x = t_rij(ipair, 0); + double y = t_rij(ipair, 1); + double z = t_rij(ipair, 2); + ipair++; + auto dplist_r=Kokkos::subview(k_dplist_r,ipair-1,Kokkos::ALL, Kokkos::ALL); + + double r = sqrt(x * x + y * y + z * z); + if (r < SMALL) continue; + + for (int ti = 0; ti < m_idxu_count; ti++) { + ulist_r[ti] = 0.0; + ulist_i[ti] = 0.0; + } + + compute_uarray_recursive(x, y, z, r, twolmax, ulist_r, ulist_i, m_idxu_block, m_rootpq); + + ///////// compute_carray_wD //////// + { + double rvec[3]; + rvec[0] = x; + rvec[1] = y; + rvec[2] = z; + + for (int i=0;i; +#ifdef LMP_KOKKOS_GPU +template class MLIAP_SO3Kokkos; +#endif +} diff --git a/src/KOKKOS/mliap_so3_kokkos.h b/src/KOKKOS/mliap_so3_kokkos.h new file mode 100644 index 0000000000..164fb3b24d --- /dev/null +++ b/src/KOKKOS/mliap_so3_kokkos.h @@ -0,0 +1,142 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) + ------------------------------------------------------------------------- */ + +#ifndef LMP_MLIAP_SO3_KOKKOS_H +#define LMP_MLIAP_SO3_KOKKOS_H + +#include "kokkos_type.h" +#include "pointers.h" + +namespace LAMMPS_NS { + +template class MLIAP_SO3Kokkos : protected Pointers { + public: + MLIAP_SO3Kokkos(LAMMPS *, double vrcut, int vlmax, int vnmax, double valpha); + MLIAP_SO3Kokkos(LAMMPS *_lmp) : Pointers(_lmp){}; + ~MLIAP_SO3Kokkos() override; + + void init(); + double memory_usage(); + + using MemoryDeviceType = typename KKDevice::value; + using float_1d = Kokkos::View; + using float_2d = Kokkos::View; + using float_3d = Kokkos::View; + using float_4d = Kokkos::View; + using int_1d = Kokkos::View; + + int ncoeff; + float_2d m_plist_r; + float_3d k_dplist_r; + + private: + double alloc_init, alloc_arrays; + int_1d m_ellpl1, m_ellm1; + float_1d m_pfac, m_Ylms; + int m_pfac_l1, m_pfac_l2; + float_1d m_dfac0, m_dfac1, m_dfac2, m_dfac3, m_dfac4, m_dfac5; + int m_dfac_l1, m_dfac_l2; + double m_rcut, m_alpha; + int m_lmax, m_nmax, m_Nmax; + float_1d m_g_array, m_w; + float_1d m_rootpq; + int_1d m_idxu_block, m_idxylm; + int m_idxu_count, m_idxy_count; + int m_numYlms; + + float_3d m_clisttot_r, m_clisttot_i; + float_1d m_rip_array, m_rip_darray; + float_1d m_sbes_array, m_sbes_darray; + int m_init_arrays; + + float_2d m_ulist_r, m_ulist_i; + + float_3d m_dYlm_r, m_dYlm_i; + float_4d m_dclist; + + // arguments for operators. + int_1d t_numneighs; + int_1d t_jelems; + float_1d t_wjelem; + float_2d t_rij; + int t_nmax; + int t_lmax; + double t_rcut; + double t_alpha; + int_1d t_ij; + + static constexpr bigint m_temp_memory_size = 512 * 1024 * 1024; + int m_chunk_size; + + public: + void spectrum(int nlocal, DAT::tdual_int_1d numneighs, DAT::tdual_int_1d jelems, + DAT::tdual_float_1d wjelem, DAT::tdual_float_2d rij, DAT::tdual_int_1d k_ij, + int nmax, int lmax, double rcut, double alpha, int totaln, int ncoefs); + struct MLIAPSO3SpectrumTag {}; + KOKKOS_FUNCTION + void operator()(const MLIAPSO3SpectrumTag &, int ii) const; + + void spectrum_dxdr(int nlocal, DAT::tdual_int_1d numneighs, DAT::tdual_int_1d jelems, + DAT::tdual_float_1d wjelem, DAT::tdual_float_2d rij, DAT::tdual_int_1d k_ij, + int nmax, int lmax, double rcut, double alpha, bigint npairs, int ncoefs); + struct MLIAPSO3SpectrumDXDRTag {}; + KOKKOS_FUNCTION + void operator()(const MLIAPSO3SpectrumDXDRTag &, int ii) const; + + KOKKOS_FUNCTION + double Cosine(double Rij, double Rc) const; + KOKKOS_FUNCTION + double CosinePrime(double Rij, double Rc) const; + KOKKOS_FUNCTION + double compute_sfac(double r, double rcut) const; + KOKKOS_FUNCTION + double compute_dsfac(double r, double rcut) const; + + struct MLIAPSO3GetSBESArrayTag {}; + KOKKOS_FUNCTION + void operator()(const MLIAPSO3GetSBESArrayTag &, int ii) const; + + struct MLIAPSO3GetRipArrayTag {}; + KOKKOS_FUNCTION + void operator()(const MLIAPSO3GetRipArrayTag &, int ii) const; + + void init_arrays(int nlocal, int ncoefs); + void init_garray(int nmax, int lmax, double rcut, double alpha, double *w, int lw1, + double *g_array, int lg2); + + template + KOKKOS_FUNCTION void compute_uarray_recursive(double x, double y, double z, double r, int twol, + UlistView ulist_r, UlistView ulist_i, + int_1d idxu_block, float_1d rootpqarray) const; + void compute_ncoeff(); + + KOKKOS_FUNCTION + static int get_sum(int istart, int iend, int id, int imult); + + double compute_g(double r, int n, int nmax, double rcut, double *w, int lw1); + double phi(double r, int alpha, double rcut); + + template + KOKKOS_FUNCTION void compute_pi(int nmax, int lmax, ViewType clisttot_r, ViewType clisttot_i, + int lcl2, float_2d plist_r, int indpl) const; + + void compute_W(int nmax, double *arr); +}; + +} // namespace LAMMPS_NS + +#endif diff --git a/src/KOKKOS/pair_mliap_kokkos.cpp b/src/KOKKOS/pair_mliap_kokkos.cpp new file mode 100644 index 0000000000..5f00c7473a --- /dev/null +++ b/src/KOKKOS/pair_mliap_kokkos.cpp @@ -0,0 +1,288 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) +------------------------------------------------------------------------- */ + +#include "pair_mliap_kokkos.h" +#include "memory_kokkos.h" +#include "atom_kokkos.h" +#include "atom_masks.h" +#include "mliap_data_kokkos.h" +#include "mliap_descriptor_so3_kokkos.h" +#include "mliap_model_linear_kokkos.h" +#include "error.h" +#include "neigh_request.h" +#include "lammps.h" +#include "kokkos.h" +#include "pointers.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +template +PairMLIAPKokkos::PairMLIAPKokkos(class LAMMPS* l) : PairMLIAP(l) +{ + kokkosable = 1; + execution_space = ExecutionSpaceFromDevice::space; + datamask_modify = 0; +} + +/* ---------------------------------------------------------------------- */ + +template +PairMLIAPKokkos::~PairMLIAPKokkos() +{ + memoryKK->destroy_kokkos(k_map, map); + memoryKK->destroy_kokkos(k_cutsq, cutsq); + memoryKK->destroy_kokkos(k_setflag, setflag); + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->destroy_kokkos(k_vatom,vatom); + allocated = 0; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairMLIAPKokkos::compute(int eflag, int vflag) +{ + atomKK->sync(Host,F_MASK | ENERGY_MASK | VIRIAL_MASK); + atomKK->sync(execution_space,X_MASK | TYPE_MASK ); + MLIAPDataKokkos *k_data = (MLIAPDataKokkos*)(data); + + int is_kokkos_model = (dynamic_cast*>(model)) != nullptr; + int is_kokkos_descriptor = (dynamic_cast*>(descriptor)) != nullptr; + auto model_space = is_kokkos_model ? execution_space : Host; + auto descriptor_space = is_kokkos_descriptor? execution_space : Host; + + // consistency checks + if (data->ndescriptors != model->ndescriptors) + error->all(FLERR, "Incompatible model and descriptor descriptor count"); + + if (data->nelements != model->nelements) + error->all(FLERR, "Incompatible model and descriptor element count"); + + ev_init(eflag, vflag); + if (eflag_atom && (int)k_eatom.h_view.extent(0) < maxeatom) { + memoryKK->destroy_kokkos(k_eatom,eatom); + memoryKK->create_kokkos(k_eatom,eatom,maxeatom,"pair:eatom"); + } + + if (vflag_atom && (int)k_vatom.h_view.extent(0) < maxeatom) { + memoryKK->destroy_kokkos(k_vatom,vatom); + memoryKK->create_kokkos(k_vatom,vatom,maxeatom,6,"pair:eatom"); + } + + data->generate_neighdata(list, eflag, vflag); + + // compute descriptors, if needed + if (model->nonlinearflag || eflag) { + k_data->sync(descriptor_space, NUMNEIGHS_MASK | IATOMS_MASK | IELEMS_MASK | JATOMS_MASK | JELEMS_MASK | RIJ_MASK ); + descriptor->compute_descriptors(data); + if (!is_kokkos_descriptor) + k_data->modified(descriptor_space, DESCRIPTORS_MASK); + } + + // compute E_i and beta_i = dE_i/dB_i for all i in list + k_data->sync(model_space, IELEMS_MASK | DESCRIPTORS_MASK); + model->compute_gradients(data); + k_data->modified(model_space, BETAS_MASK); + if (eflag_atom) + k_data->modified(model_space, EATOMS_MASK); + e_tally(data); + + // calculate force contributions beta_i*dB_i/dR_j + k_data->sync(descriptor_space, NUMNEIGHS_MASK | IATOMS_MASK | IELEMS_MASK | BETAS_MASK | JATOMS_MASK | JELEMS_MASK | RIJ_MASK ); + descriptor->compute_forces(data); + + if (evflag) { + atomKK->modified(descriptor_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); + atomKK->sync(execution_space,F_MASK | ENERGY_MASK | VIRIAL_MASK); + } else { + atomKK->modified(descriptor_space,F_MASK); + atomKK->sync(execution_space,F_MASK); + } + + // calculate stress + if (vflag_fdotr) { + pair_virial_fdotr_compute(this); + } +} + +/* ---------------------------------------------------------------------- */ + +template +void PairMLIAPKokkos::allocate() +{ + int n = atom->ntypes; + memoryKK->destroy_kokkos(k_map, map); + memoryKK->destroy_kokkos(k_cutsq, cutsq); + memoryKK->destroy_kokkos(k_setflag, setflag); + memoryKK->create_kokkos(k_map, map, n+1, "pair_mliap:map"); + memoryKK->create_kokkos(k_cutsq, cutsq, n+1, n+1, "pair_mliap:cutsq"); + memoryKK->create_kokkos(k_setflag, setflag, n+1, n+1, "pair_mliap:setflag"); + + // this is for the base class so it doesn't double delete + allocated = 1; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairMLIAPKokkos::settings(int narg, char ** arg) +{ + PairMLIAP::settings(narg, arg); + int iarg=0; + while (iarg < narg) { + if (strcmp(arg[iarg],"model") == 0) { + if (strcmp(arg[iarg+1],"linear") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + delete model; + model = new MLIAPModelLinearKokkos(lmp,arg[iarg+2]); + iarg += 3; + } else + iarg += 2; + } else if (strcmp(arg[iarg],"descriptor") == 0) { + if (strcmp(arg[iarg+1],"so3") == 0) { + if (iarg+3 > narg) error->all(FLERR,"Illegal pair_style mliap command"); + delete descriptor; + descriptor = new MLIAPDescriptorSO3Kokkos(lmp,arg[iarg+2]); + iarg += 3; + } else + iarg ++; + } else + iarg++; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void PairMLIAPKokkos::coeff(int narg, char **arg) { + if (narg < 3) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) { + PairMLIAP::allocate(); + allocate(); + } + + char* type1 = arg[0]; + char* type2 = arg[1]; + char** elemtypes = &arg[2]; + + // insure I,J args are * * + + if (strcmp(type1,"*") != 0 || strcmp(type2,"*") != 0) + error->all(FLERR,"Incorrect args for pair coefficients"); + + // read args that map atom types to elements + // map[i] = which element the Ith atom type is, -1 if not mapped + // map[0] is not used + + for (int i = 1; i <= atom->ntypes; i++) { + char* elemname = elemtypes[i-1]; + int jelem; + for (jelem = 0; jelem < descriptor->nelements; jelem++) + if (strcmp(elemname,descriptor->elements[jelem]) == 0) + break; + + if (jelem < descriptor->nelements) + map[i] = jelem; + else if (strcmp(elemname,"NULL") == 0) map[i] = -1; + else error->all(FLERR,"Incorrect args for pair coefficients"); + } + k_map.modify(); + k_map.sync(); + + auto h_cutsq=k_cutsq.template view(); + for (int itype=1; itype <= atom->ntypes; ++itype) + for (int jtype=1; jtype <= atom->ntypes; ++jtype) + h_cutsq(itype,jtype) = descriptor->cutsq[map[itype]][map[jtype]]; + k_cutsq.modify(); + k_cutsq.sync(); + + // clear setflag since coeff() called once with I,J = * * + + int n = atom->ntypes; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + // set setflag i,j for type pairs where both are mapped to elements + + int count = 0; + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + if (map[i] >= 0 && map[j] >= 0) { + setflag[i][j] = 1; + count++; + } + k_setflag.modify(); + k_setflag.sync(); + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); + + // set up model, descriptor, and mliap data structures + model->init(); + descriptor->init(); + int gradgradflag = -1; + delete data; + data = new MLIAPDataKokkos(lmp, gradgradflag, map, model, descriptor, this); + data->init(); +} + +/* ---------------------------------------------------------------------- + add energies to eng_vdwl and per-atom energy +------------------------------------------------------------------------- */ + +template +void PairMLIAPKokkos::e_tally(MLIAPData* data) +{ + if (eflag_global) eng_vdwl += data->energy; + if (eflag_atom) { + MLIAPDataKokkos *k_data = static_cast*>(data); + k_data->sync(execution_space, IATOMS_MASK | EATOMS_MASK); + auto d_iatoms = k_data->k_iatoms.template view(); + auto d_eatoms = k_data->k_eatoms.template view(); + auto d_eatom = k_eatom.template view(); + Kokkos::parallel_for(data->nlistatoms, KOKKOS_LAMBDA (int ii) { + d_eatom(d_iatoms(ii)) = d_eatoms(ii); + }); + k_eatom.modify(); + // This sync has to be here for the hybrid pair type + k_eatom.sync(); + } +} + +/* ---------------------------------------------------------------------- */ + +template +void PairMLIAPKokkos::init_style() +{ + + PairMLIAP::init_style(); + auto request = neighbor->find_request(this); + request->set_kokkos_host(std::is_same::value && + !std::is_same::value); + request->set_kokkos_device(std::is_same::value); +} + +/* ---------------------------------------------------------------------- */ + +namespace LAMMPS_NS { +template class PairMLIAPKokkos; +#ifdef LMP_KOKKOS_GPU +template class PairMLIAPKokkos; +#endif +} diff --git a/src/KOKKOS/pair_mliap_kokkos.h b/src/KOKKOS/pair_mliap_kokkos.h new file mode 100644 index 0000000000..1acac6bf24 --- /dev/null +++ b/src/KOKKOS/pair_mliap_kokkos.h @@ -0,0 +1,67 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS Development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Matt Bettencourt (NVIDIA) +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS +// clang-format off +PairStyle(mliap/kk,PairMLIAPKokkos); +PairStyle(mliap/kk/device,PairMLIAPKokkos); +PairStyle(mliap/kk/host,PairMLIAPKokkos); +// clang-format off +#else + +#ifndef LMP_PAIR_MLIAP_KOKKOS_H +#define LMP_PAIR_MLIAP_KOKKOS_H + +#include "pair_mliap.h" +#include "pair_kokkos.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +template +class PairMLIAPKokkos : public PairMLIAP { +public: + typedef DeviceType device_type; + typedef ArrayTypes AT; + + PairMLIAPKokkos(class LAMMPS*); + ~PairMLIAPKokkos(); + void settings(int narg, char ** arg); + void init_style(); + + void compute(int, int); + void e_tally(MLIAPData* data); + + void allocate(); + + void coeff(int narg, char **arg); + typename AT::t_x_array_randomread x; + typename AT::t_x_array_randomread v; + typename AT::t_f_array f; + DAT::tdual_int_1d k_map; + DAT::tdual_double_2d k_cutsq; + DAT::tdual_int_2d k_setflag; + DAT::tdual_efloat_1d k_eatom; + DAT::tdual_double_2d k_vatom; + + + friend void pair_virial_fdotr_compute(PairMLIAPKokkos*); +}; + +} +#endif +#endif diff --git a/src/ML-IAP/mliap_data.h b/src/ML-IAP/mliap_data.h index 28c50a440c..562eb3e850 100644 --- a/src/ML-IAP/mliap_data.h +++ b/src/ML-IAP/mliap_data.h @@ -26,8 +26,8 @@ class MLIAPData : protected Pointers { ~MLIAPData() override; void init(); - void generate_neighdata(class NeighList *, int = 0, int = 0); - void grow_neigharrays(); + virtual void generate_neighdata(class NeighList *, int = 0, int = 0); + virtual void grow_neigharrays(); double memory_usage(); int size_array_rows, size_array_cols; @@ -77,7 +77,7 @@ class MLIAPData : protected Pointers { int vflag; // indicates if virial is needed class PairMLIAP *pairmliap; // access to pair tally functions - private: + protected: class MLIAPModel *model; class MLIAPDescriptor *descriptor; diff --git a/src/ML-IAP/mliap_descriptor.h b/src/ML-IAP/mliap_descriptor.h index 9775e55b5a..7f05c6804a 100644 --- a/src/ML-IAP/mliap_descriptor.h +++ b/src/ML-IAP/mliap_descriptor.h @@ -18,7 +18,7 @@ namespace LAMMPS_NS { -class MLIAPDescriptor : protected Pointers { +class MLIAPDescriptor : virtual protected Pointers { public: MLIAPDescriptor(LAMMPS *); ~MLIAPDescriptor() override; diff --git a/src/ML-IAP/mliap_descriptor_snap.cpp b/src/ML-IAP/mliap_descriptor_snap.cpp index 7b521b065d..b8f9aa087b 100644 --- a/src/ML-IAP/mliap_descriptor_snap.cpp +++ b/src/ML-IAP/mliap_descriptor_snap.cpp @@ -36,7 +36,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *_lmp, char *paramfilename) : MLIAPDescriptor(_lmp) +MLIAPDescriptorSNAP::MLIAPDescriptorSNAP(LAMMPS *_lmp, char *paramfilename) : Pointers(_lmp), MLIAPDescriptor(_lmp) { radelem = nullptr; wjelem = nullptr; diff --git a/src/ML-IAP/mliap_descriptor_so3.cpp b/src/ML-IAP/mliap_descriptor_so3.cpp index dbf383917f..96a589eadd 100644 --- a/src/ML-IAP/mliap_descriptor_so3.cpp +++ b/src/ML-IAP/mliap_descriptor_so3.cpp @@ -10,6 +10,7 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ + /* ---------------------------------------------------------------------- Contributing authors: Byungkyun Kang (University of Nevada, Las Vegas) ------------------------------------------------------------------------- */ @@ -34,7 +35,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -MLIAPDescriptorSO3::MLIAPDescriptorSO3(LAMMPS *lmp, char *paramfilename) : MLIAPDescriptor(lmp) +MLIAPDescriptorSO3::MLIAPDescriptorSO3(LAMMPS *lmp, char *paramfilename) : Pointers(lmp), MLIAPDescriptor(lmp) { radelem = nullptr; wjelem = nullptr; diff --git a/src/ML-IAP/mliap_descriptor_so3.h b/src/ML-IAP/mliap_descriptor_so3.h index 8f55374888..94785f63e5 100644 --- a/src/ML-IAP/mliap_descriptor_so3.h +++ b/src/ML-IAP/mliap_descriptor_so3.h @@ -18,7 +18,7 @@ namespace LAMMPS_NS { -class MLIAPDescriptorSO3 : public MLIAPDescriptor { +class MLIAPDescriptorSO3 : public MLIAPDescriptor, virtual protected Pointers { public: MLIAPDescriptorSO3(LAMMPS *, char *); diff --git a/src/ML-IAP/mliap_model_linear.cpp b/src/ML-IAP/mliap_model_linear.cpp index e966afd4af..4d8500ae43 100644 --- a/src/ML-IAP/mliap_model_linear.cpp +++ b/src/ML-IAP/mliap_model_linear.cpp @@ -57,22 +57,17 @@ void MLIAPModelLinear::compute_gradients(MLIAPData* data) for (int ii = 0; ii < data->nlistatoms; ii++) { const int ielem = data->ielems[ii]; - double* coeffi = coeffelem[ielem]; + double const* coeffi = coeffelem[ielem]; for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) data->betas[ii][icoeff] = coeffi[icoeff+1]; // add in contributions to global and per-atom energy // this is optional and has no effect on force calculation - if (data->eflag) { - // energy of atom I - - double* coeffi = coeffelem[ielem]; double etmp = coeffi[0]; // E_i = beta.B_i - for (int icoeff = 0; icoeff < data->ndescriptors; icoeff++) etmp += coeffi[icoeff+1]*data->descriptors[ii][icoeff]; diff --git a/src/ML-IAP/mliap_so3.cpp b/src/ML-IAP/mliap_so3.cpp index 7b9a683253..916b8cfcf5 100644 --- a/src/ML-IAP/mliap_so3.cpp +++ b/src/ML-IAP/mliap_so3.cpp @@ -1,18 +1,18 @@ /* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/, Sandia National Laboratories - LAMMPS development team: developers@lammps.org + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. - See the README file in the top-level LAMMPS directory. + See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: Byungkyun Kang (University of Nevada, Las Vegas) + Contributing authors: Byungkyun Kang (University of Nevada, Las Vegas) ------------------------------------------------------------------------- */ #include "mliap_so3.h" @@ -84,7 +84,6 @@ MLIAP_SO3::MLIAP_SO3(LAMMPS *lmp, double vrcut, int vlmax, int vnmax, double val m_dYlm_i = nullptr; m_dplist_r = nullptr; - m_dplist_i = nullptr; m_dclist_r = nullptr; m_dclist_i = nullptr; @@ -145,7 +144,6 @@ MLIAP_SO3::~MLIAP_SO3() memory->destroy(m_dYlm_i); memory->destroy(m_dplist_r); - memory->destroy(m_dplist_i); memory->destroy(m_dclist_r); memory->destroy(m_dclist_i); @@ -846,8 +844,6 @@ void MLIAP_SO3::spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem totali = totaln * ncoefs * 3; memory->destroy(m_dplist_r); memory->create(m_dplist_r, totali, "MLIAP_SO3:m_dplist_r"); - memory->destroy(m_dplist_i); - memory->create(m_dplist_i, totali, "MLIAP_SO3:m_dplist_i"); alloc_arrays += 2.0 * totali * sizeof(double); get_sbes_array(nlocal, numneighs, rij, lmax, rcut, alpha); @@ -977,14 +973,11 @@ void MLIAP_SO3::spectrum_dxdr(int nlocal, int *numneighs, int *jelems, double *w totali = totaln * ncoefs * 3; memory->destroy(m_dplist_r); memory->create(m_dplist_r, totali, "MLIAP_SO3:m_dplist_r"); - memory->destroy(m_dplist_i); - memory->create(m_dplist_i, totali, "MLIAP_SO3:m_dplist_i"); totali = npairs * ncoefs * 3; for (int i = 0; i < totali; i++) { m_dplist_r[i] = 0.0; - m_dplist_i[i] = 0.0; } numps = nmax * (nmax + 1) * (lmax + 1) / 2; diff --git a/src/ML-IAP/mliap_so3.h b/src/ML-IAP/mliap_so3.h index 27d4f406e4..caf6adff76 100644 --- a/src/ML-IAP/mliap_so3.h +++ b/src/ML-IAP/mliap_so3.h @@ -14,7 +14,6 @@ class MLIAP_SO3 : protected Pointers { ~MLIAP_SO3() override; void init(); - double memory_usage(); int ncoeff; @@ -46,7 +45,7 @@ class MLIAP_SO3 : protected Pointers { double *m_ulist_r, *m_ulist_i; double *m_Ylms_r, *m_Ylms_i, *m_dYlm_r, *m_dYlm_i; - double *m_dplist_i, *m_dclist_r, *m_dclist_i, *m_tempdp_r; + double *m_dclist_r, *m_dclist_i, *m_tempdp_r; public: void spectrum(int nlocal, int *numneighs, int *jelems, double *wjelem, double **rij, int nmax, diff --git a/src/ML-IAP/mliap_so3_math.h b/src/ML-IAP/mliap_so3_math.h index 1e26b4d223..21f7261496 100644 --- a/src/ML-IAP/mliap_so3_math.h +++ b/src/ML-IAP/mliap_so3_math.h @@ -9,6 +9,7 @@ #define LMP_MLIAP_SO3_MATH_H #include "math_eigen_impl.h" +#include namespace SO3Math { @@ -22,7 +23,7 @@ void LUPSolve(int n, double *A, double *B, int *P); using namespace MathEigen; typedef Jacobi Jacobi_v2; -int SO3Math::jacobin(int n, double const *const *mat, double *eval, double **evec) +inline int SO3Math::jacobin(int n, double const *const *mat, double *eval, double **evec) { int *midx = new int[n]; double **M = new double *[n]; @@ -48,7 +49,7 @@ int SO3Math::jacobin(int n, double const *const *mat, double *eval, double **eve return ierror; } -int SO3Math::invert_matrix(int n, double *A, double *Ainv) +inline int SO3Math::invert_matrix(int n, double *A, double *Ainv) { int i, j; @@ -85,7 +86,7 @@ int SO3Math::invert_matrix(int n, double *A, double *Ainv) return rv; } -int SO3Math::LUPdecompose(int n, double dtol, double *A, int *P) +inline int SO3Math::LUPdecompose(int n, double dtol, double *A, int *P) { int i, j, k, maxi; double maxA, Atemp; @@ -148,7 +149,7 @@ int SO3Math::LUPdecompose(int n, double dtol, double *A, int *P) return 0; } -void SO3Math::LUPSolve(int n, double *A, double *B, int *P) +inline void SO3Math::LUPSolve(int n, double *A, double *B, int *P) { int i, j; double dtemp; diff --git a/src/ML-IAP/mliap_unified.cpp b/src/ML-IAP/mliap_unified.cpp index 06cea7693f..a3727976de 100644 --- a/src/ML-IAP/mliap_unified.cpp +++ b/src/ML-IAP/mliap_unified.cpp @@ -33,7 +33,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -MLIAPDummyDescriptor::MLIAPDummyDescriptor(LAMMPS *_lmp) : MLIAPDescriptor(_lmp) {} +MLIAPDummyDescriptor::MLIAPDummyDescriptor(LAMMPS *_lmp) : Pointers(_lmp), MLIAPDescriptor(_lmp) {} MLIAPDummyDescriptor::~MLIAPDummyDescriptor() { diff --git a/src/ML-IAP/mliap_unified.h b/src/ML-IAP/mliap_unified.h index 46eeee8edd..1ee8dfd60d 100644 --- a/src/ML-IAP/mliap_unified.h +++ b/src/ML-IAP/mliap_unified.h @@ -22,7 +22,7 @@ namespace LAMMPS_NS { -class MLIAPDummyDescriptor : public MLIAPDescriptor { +class MLIAPDummyDescriptor : public MLIAPDescriptor, virtual protected Pointers { public: MLIAPDummyDescriptor(LAMMPS *); ~MLIAPDummyDescriptor() override; diff --git a/src/ML-IAP/pair_mliap.cpp b/src/ML-IAP/pair_mliap.cpp index 2a7757b53d..38d98510ce 100644 --- a/src/ML-IAP/pair_mliap.cpp +++ b/src/ML-IAP/pair_mliap.cpp @@ -135,10 +135,10 @@ void PairMLIAP::settings(int narg, char ** arg) // set flags for required keywords - int modelflag = 0; - int descriptorflag = 0; delete model; + model = nullptr; delete descriptor; + descriptor = nullptr; // process keywords @@ -147,7 +147,7 @@ void PairMLIAP::settings(int narg, char ** arg) while (iarg < narg) { if (strcmp(arg[iarg],"model") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model", error); - if (modelflag) error->all(FLERR,"Illegal multiple pair_style mliap model definition"); + if (model != nullptr) error->all(FLERR,"Illegal multiple pair_style mliap model definition"); if (strcmp(arg[iarg+1],"linear") == 0) { if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap model linear", error); model = new MLIAPModelLinear(lmp,arg[iarg+2]); @@ -169,10 +169,9 @@ void PairMLIAP::settings(int narg, char ** arg) error->all(FLERR,"Using pair_style mliap model mliappy requires ML-IAP with python support"); #endif } else error->all(FLERR,"Unknown pair_style mliap model keyword: {}", arg[iarg]); - modelflag = 1; } else if (strcmp(arg[iarg],"descriptor") == 0) { if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor", error); - if (descriptorflag) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definition"); + if (descriptor != nullptr) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definition"); if (strcmp(arg[iarg+1],"sna") == 0) { if (iarg+3 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap descriptor sna", error); descriptor = new MLIAPDescriptorSNAP(lmp,arg[iarg+2]); @@ -183,11 +182,10 @@ void PairMLIAP::settings(int narg, char ** arg) iarg += 3; } else error->all(FLERR,"Illegal pair_style mliap command"); - descriptorflag = 1; } else if (strcmp(arg[iarg], "unified") == 0) { #ifdef MLIAP_PYTHON - if (modelflag) error->all(FLERR,"Illegal multiple pair_style mliap model definitions"); - if (descriptorflag) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definitions"); + if (model != nullptr) error->all(FLERR,"Illegal multiple pair_style mliap model definitions"); + if (descriptor != nullptr) error->all(FLERR,"Illegal multiple pair_style mliap descriptor definitions"); if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "pair_style mliap unified", error); MLIAPBuildUnified_t build = build_unified(arg[iarg+1], data, lmp); if (iarg+3 > narg) { @@ -199,8 +197,6 @@ void PairMLIAP::settings(int narg, char ** arg) iarg += 3; model = build.model; descriptor = build.descriptor; - modelflag = 1; - descriptorflag = 1; #else error->all(FLERR,"Using pair_style mliap unified requires ML-IAP with python support"); #endif @@ -208,7 +204,7 @@ void PairMLIAP::settings(int narg, char ** arg) error->all(FLERR,"Unknown pair_style mliap keyword: {}", arg[iarg]); } - if (modelflag == 0 || descriptorflag == 0) + if (model == nullptr || descriptor == nullptr) error->all(FLERR,"Incomplete pair_style mliap setup: need model and descriptor, or unified"); }