From 460e90ac5f75e69ccb04cd76b5d9f7a1ffca89b1 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 13 Jan 2021 16:48:27 -0500 Subject: [PATCH 1/6] Move Cabana interface to Kokkos --- .../{CabanaMD/ElementCabana.h => Kokkos/ElementKokkos.h} | 0 .../ElementCabana_impl.h => Kokkos/ElementKokkos_impl.h} | 0 src/libnnpif/{CabanaMD/ModeCabana.h => Kokkos/ModeKokkos.h} | 2 +- .../{CabanaMD/ModeCabana_impl.h => Kokkos/ModeKokkos_impl.h} | 0 src/libnnpif/{CabanaMD/typesCabana.h => Kokkos/typesKokkos.h} | 0 5 files changed, 1 insertion(+), 1 deletion(-) rename src/libnnpif/{CabanaMD/ElementCabana.h => Kokkos/ElementKokkos.h} (100%) rename src/libnnpif/{CabanaMD/ElementCabana_impl.h => Kokkos/ElementKokkos_impl.h} (100%) rename src/libnnpif/{CabanaMD/ModeCabana.h => Kokkos/ModeKokkos.h} (99%) rename src/libnnpif/{CabanaMD/ModeCabana_impl.h => Kokkos/ModeKokkos_impl.h} (100%) rename src/libnnpif/{CabanaMD/typesCabana.h => Kokkos/typesKokkos.h} (100%) diff --git a/src/libnnpif/CabanaMD/ElementCabana.h b/src/libnnpif/Kokkos/ElementKokkos.h similarity index 100% rename from src/libnnpif/CabanaMD/ElementCabana.h rename to src/libnnpif/Kokkos/ElementKokkos.h diff --git a/src/libnnpif/CabanaMD/ElementCabana_impl.h b/src/libnnpif/Kokkos/ElementKokkos_impl.h similarity index 100% rename from src/libnnpif/CabanaMD/ElementCabana_impl.h rename to src/libnnpif/Kokkos/ElementKokkos_impl.h diff --git a/src/libnnpif/CabanaMD/ModeCabana.h b/src/libnnpif/Kokkos/ModeKokkos.h similarity index 99% rename from src/libnnpif/CabanaMD/ModeCabana.h rename to src/libnnpif/Kokkos/ModeKokkos.h index 2374c8e01..f64d784be 100644 --- a/src/libnnpif/CabanaMD/ModeCabana.h +++ b/src/libnnpif/Kokkos/ModeKokkos.h @@ -18,7 +18,7 @@ #ifndef MODE_CABANA_H #define MODE_CABANA_H -#include "./ElementCabana.h" +#include "ElementCabana.h" #include "typesCabana.h" #include diff --git a/src/libnnpif/CabanaMD/ModeCabana_impl.h b/src/libnnpif/Kokkos/ModeKokkos_impl.h similarity index 100% rename from src/libnnpif/CabanaMD/ModeCabana_impl.h rename to src/libnnpif/Kokkos/ModeKokkos_impl.h diff --git a/src/libnnpif/CabanaMD/typesCabana.h b/src/libnnpif/Kokkos/typesKokkos.h similarity index 100% rename from src/libnnpif/CabanaMD/typesCabana.h rename to src/libnnpif/Kokkos/typesKokkos.h From ee87fad41662e82eee1ae9f4b6954265c48a7776 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 13 Jan 2021 16:50:25 -0500 Subject: [PATCH 2/6] Rename everything to Kokkos; remove Cabana features from main compute kernels --- src/libnnpif/Kokkos/ElementKokkos.h | 24 ++++----- src/libnnpif/Kokkos/ElementKokkos_impl.h | 30 +++++------ src/libnnpif/Kokkos/ModeKokkos.h | 62 ++++++++++----------- src/libnnpif/Kokkos/ModeKokkos_impl.h | 68 ++++++++++++------------ src/libnnpif/Kokkos/typesKokkos.h | 4 +- 5 files changed, 89 insertions(+), 99 deletions(-) diff --git a/src/libnnpif/Kokkos/ElementKokkos.h b/src/libnnpif/Kokkos/ElementKokkos.h index 69328dd8d..09180dece 100644 --- a/src/libnnpif/Kokkos/ElementKokkos.h +++ b/src/libnnpif/Kokkos/ElementKokkos.h @@ -15,10 +15,10 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#ifndef ELEMENT_CABANA_H -#define ELEMENT_CABANA_H +#ifndef ELEMENT_KOKKOS_H +#define ELEMENT_KOKKOS_H -#include "typesCabana.h" +#include "typesKokkos.h" #include "CutoffFunction.h" #include "Element.h" @@ -31,16 +31,16 @@ namespace nnp { -/// Derived Cabana class for element-specific data. -class ElementCabana : public Element +/// Derived Kokkos class for element-specific data. +class ElementKokkos : public Element { public: /** Constructor using index. */ - ElementCabana( std::size_t const index ); + ElementKokkos( std::size_t const index ); /** Destructor. */ - ~ElementCabana(); + ~ElementKokkos(); /** Add one symmetry function. * @@ -275,7 +275,7 @@ class ElementCabana : public Element template inline std::size_t -ElementCabana::numSymmetryFunctions( int attype, +ElementKokkos::numSymmetryFunctions( int attype, h_t_int h_numSFperElem ) const { return h_numSFperElem( attype ); @@ -283,7 +283,7 @@ ElementCabana::numSymmetryFunctions( int attype, template inline void -ElementCabana::setScalingType( ScalingType scalingType, +ElementKokkos::setScalingType( ScalingType scalingType, std::string statisticsLine, double Smin, double Smax, t_SFscaling SFscaling, int attype, int k ) const @@ -324,7 +324,7 @@ ElementCabana::setScalingType( ScalingType scalingType, template inline std::string -ElementCabana::scalingLine( ScalingType scalingType, +ElementKokkos::scalingLine( ScalingType scalingType, t_SFscaling SFscaling, int attype, int k ) const { @@ -337,7 +337,7 @@ ElementCabana::scalingLine( ScalingType scalingType, } template -inline double ElementCabana::unscale( int attype, double value, int k, +inline double ElementKokkos::unscale( int attype, double value, int k, t_SFscaling SFscaling ) { double scalingType = SFscaling( attype, k, 7 ); @@ -377,6 +377,6 @@ inline double ElementCabana::unscale( int attype, double value, int k, } -#include "ElementCabana_impl.h" +#include "ElementKokkos_impl.h" #endif diff --git a/src/libnnpif/Kokkos/ElementKokkos_impl.h b/src/libnnpif/Kokkos/ElementKokkos_impl.h index 3feae0507..f952b608a 100644 --- a/src/libnnpif/Kokkos/ElementKokkos_impl.h +++ b/src/libnnpif/Kokkos/ElementKokkos_impl.h @@ -28,7 +28,7 @@ using namespace std; namespace nnp { -ElementCabana::ElementCabana( size_t const _index ) +ElementKokkos::ElementKokkos( size_t const _index ) : Element() { index = _index; @@ -36,10 +36,10 @@ ElementCabana::ElementCabana( size_t const _index ) neuralNetwork = NULL; } -ElementCabana::~ElementCabana() {} +ElementKokkos::~ElementKokkos() {} template -void ElementCabana::addSymmetryFunction( string const ¶meters, +void ElementKokkos::addSymmetryFunction( string const ¶meters, vector elementStrings, int attype, t_SF SF, double convLength, h_t_int h_numSFperElem ) @@ -166,7 +166,7 @@ void ElementCabana::addSymmetryFunction( string const ¶meters, } template -void ElementCabana::sortSymmetryFunctions( t_SF SF, h_t_int h_numSFperElem, +void ElementKokkos::sortSymmetryFunctions( t_SF SF, h_t_int h_numSFperElem, int attype ) { int size = h_numSFperElem( attype ); @@ -202,7 +202,7 @@ void ElementCabana::sortSymmetryFunctions( t_SF SF, h_t_int h_numSFperElem, } template -bool ElementCabana::compareSF( t_SF SF, int attype, int index1, int index2 ) +bool ElementKokkos::compareSF( t_SF SF, int attype, int index1, int index2 ) { if ( SF( attype, index2, 0 ) < SF( attype, index1, 0 ) ) return true; // ec @@ -265,7 +265,7 @@ bool ElementCabana::compareSF( t_SF SF, int attype, int index1, int index2 ) template vector -ElementCabana::infoSymmetryFunctionParameters( t_SF SF, int attype, +ElementKokkos::infoSymmetryFunctionParameters( t_SF SF, int attype, h_t_int h_numSFperElem ) const { vector v; @@ -290,7 +290,7 @@ ElementCabana::infoSymmetryFunctionParameters( t_SF SF, int attype, template vector -ElementCabana::infoSymmetryFunctionScaling( ScalingType scalingType, t_SF SF, +ElementKokkos::infoSymmetryFunctionScaling( ScalingType scalingType, t_SF SF, t_SFscaling SFscaling, int attype, h_t_int h_numSFperElem ) const { @@ -305,7 +305,7 @@ ElementCabana::infoSymmetryFunctionScaling( ScalingType scalingType, t_SF SF, } template -void ElementCabana::setupSymmetryFunctionGroups( t_SF SF, +void ElementKokkos::setupSymmetryFunctionGroups( t_SF SF, t_SFGmemberlist SFGmemberlist, int attype, h_t_int h_numSFperElem, h_t_int h_numSFGperElem, @@ -387,7 +387,7 @@ void ElementCabana::setupSymmetryFunctionGroups( t_SF SF, template vector -ElementCabana::infoSymmetryFunctionGroups( t_SF SF, t_SFGmemberlist SFGmemberlist, +ElementKokkos::infoSymmetryFunctionGroups( t_SF SF, t_SFGmemberlist SFGmemberlist, int attype, h_t_int h_numSFGperElem ) const { vector v; @@ -409,7 +409,7 @@ ElementCabana::infoSymmetryFunctionGroups( t_SF SF, t_SFGmemberlist SFGmemberlis } template -void ElementCabana::setCutoffFunction( CutoffFunction::CutoffType const cutoffType, +void ElementKokkos::setCutoffFunction( CutoffFunction::CutoffType const cutoffType, double const cutoffAlpha, t_SF SF, int attype, h_t_int h_numSFperElem ) { @@ -422,7 +422,7 @@ void ElementCabana::setCutoffFunction( CutoffFunction::CutoffType const cutoffTy } template -void ElementCabana::setScaling( ScalingType scalingType, +void ElementKokkos::setScaling( ScalingType scalingType, vector const &statisticsLine, double Smin, double Smax, t_SF SF, t_SFscaling SFscaling, int attype, h_t_int h_numSFperElem ) const @@ -442,7 +442,7 @@ void ElementCabana::setScaling( ScalingType scalingType, } template -size_t ElementCabana::getMinNeighbors( int attype, t_SF SF, int nSF ) const +size_t ElementKokkos::getMinNeighbors( int attype, t_SF SF, int nSF ) const { // get max number of minNeighbors size_t global_minNeighbors = 0; @@ -462,7 +462,7 @@ size_t ElementCabana::getMinNeighbors( int attype, t_SF SF, int nSF ) const } template -double ElementCabana::getMinCutoffRadius( t_SF SF, int attype, +double ElementKokkos::getMinCutoffRadius( t_SF SF, int attype, h_t_int h_numSFperElem ) const { double minCutoffRadius = numeric_limits::max(); @@ -474,7 +474,7 @@ double ElementCabana::getMinCutoffRadius( t_SF SF, int attype, } template -double ElementCabana::getMaxCutoffRadius( t_SF SF, int attype, +double ElementKokkos::getMaxCutoffRadius( t_SF SF, int attype, h_t_int h_numSFperElem ) const { double maxCutoffRadius = 0.0; @@ -487,7 +487,7 @@ double ElementCabana::getMaxCutoffRadius( t_SF SF, int attype, // TODO:add functionality /* -void ElementCabana::updateSymmetryFunctionStatistics +void ElementKokkos::updateSymmetryFunctionStatistics */ } diff --git a/src/libnnpif/Kokkos/ModeKokkos.h b/src/libnnpif/Kokkos/ModeKokkos.h index f64d784be..33283e580 100644 --- a/src/libnnpif/Kokkos/ModeKokkos.h +++ b/src/libnnpif/Kokkos/ModeKokkos.h @@ -15,13 +15,12 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#ifndef MODE_CABANA_H -#define MODE_CABANA_H +#ifndef MODE_KOKKOS_H +#define MODE_KOKKOS_H -#include "ElementCabana.h" -#include "typesCabana.h" +#include "ElementKokkos.h" +#include "typesKokkos.h" -#include #include #include "CutoffFunction.h" @@ -36,14 +35,14 @@ namespace nnp { -/** Derived Cabana main NNP class. +/** Derived Kokkos main NNP class. * * The main n2p2 functions for computing energies and forces are replaced - * to use the Kokkos and Cabana libraries. Most setup functions are + * to use the Kokkos library. Most setup functions are * overridden; some are replaced as needed to work within device kernels. */ template -class ModeCabana : public Mode +class ModeKokkos : public Mode { public: @@ -162,32 +161,28 @@ class ModeCabana : public Mode /** Calculate forces for all atoms in given structure. * - * @param[in] x Cabana slice of atom positions. - * @param[in] f Cabana slice of atom forces. - * @param[in] type Cabana slice of atom types. - * @param[in] dEdG Cabana slice of the derivative of energy with respect + * @param[in] x Kokkos View of atom positions. + * @param[in] f Kokkos View of atom forces. + * @param[in] type Kokkos View of atom types. + * @param[in] dEdG Kokkos View of the derivative of energy with respect to symmetry functions per atom. * @param[in] N_local Number of atoms. - * @param[in] neigh_op Cabana tag for neighbor parallelism. - * @param[in] angle_op Cabana tag for angular parallelism. * * Combine intermediate results from symmetry function and neural network * computation to atomic forces. Results are stored in f. */ template + class t_slice_dEdG, class t_neigh_list> void calculateForces(t_slice_x x, t_slice_f f, t_slice_type type, - t_slice_dEdG dEdG, t_neigh_list neigh_list, int N_local, - t_neigh_parallel neigh_op, t_angle_parallel angle_op); + t_slice_dEdG dEdG, t_neigh_list neigh_list, int N_local ); /** Calculate atomic neural networks for all atoms in given structure. * - * @param[in] type Cabana slice of atom types. - * @param[in] G Cabana slice of symmetry functions per atom. - * @param[in] dEdG Cabana slice of the derivative of energy with respect + * @param[in] type Kokkos View of atom types. + * @param[in] G Kokkos View of symmetry functions per atom. + * @param[in] dEdG Kokkos View of the derivative of energy with respect to symmetry functions per atom. - * @param[in] E Cabana slice of energy per atom. + * @param[in] E Kokkos View of energy per atom. * @param[in] N_local Number of atoms. * * The atomic energy is stored in E. @@ -201,23 +196,20 @@ class ModeCabana : public Mode /** Calculate all symmetry function groups for all atoms in given * structure. * - * @param[in] x Cabana slice of atom positions. - * @param[in] type Cabana slice of atom types. - * @param[in] G Cabana slice of symmetry functions per atom. - * @param[in] neigh_list Cabana neighbor list. + * @param[in] x Kokkos View of atom positions. + * @param[in] type Kokkos View of atom types. + * @param[in] G Kokkos View of symmetry functions per atom. + * @param[in] neigh_list neighbor list. * @param[in] N_local Number of atoms. - * @param[in] neigh_op Cabana tag for neighbor parallelism. - * @param[in] angle_op Cabana tag for angular parallelism. * * Note there is no calculateSymmetryFunctions() within this derived class. * Results are stored in G. */ template + class t_neigh_list> void calculateSymmetryFunctionGroups(t_slice_x x, t_slice_type type, t_slice_G G, t_neigh_list neigh_list, - int N_local, t_neigh_parallel neigh_op, - t_angle_parallel angle_op); + int N_local); using Mode::log; @@ -272,7 +264,7 @@ class ModeCabana : public Mode ScalingType scalingType; using Mode::settings; using Mode::cutoffType; - std::vector elements; + std::vector elements; std::vector elementStrings; }; @@ -282,7 +274,7 @@ class ModeCabana : public Mode template KOKKOS_INLINE_FUNCTION void -ModeCabana::compute_cutoff(CutoffFunction::CutoffType cutoffType, +ModeKokkos::compute_cutoff(CutoffFunction::CutoffType cutoffType, double cutoffAlpha, double &fc, double &dfc, double r, double rc, bool derivative) const { @@ -313,7 +305,7 @@ ModeCabana::compute_cutoff(CutoffFunction::CutoffType cutoffType, template KOKKOS_INLINE_FUNCTION double -ModeCabana::scale(int attype, double value, int k, +ModeKokkos::scale(int attype, double value, int k, d_t_SFscaling SFscaling_) const { double scalingType = SFscaling_(attype, k, 7); @@ -342,6 +334,6 @@ ModeCabana::scale(int attype, double value, int k, } -#include "ModeCabana_impl.h" +#include "ModeKokkos_impl.h" #endif diff --git a/src/libnnpif/Kokkos/ModeKokkos_impl.h b/src/libnnpif/Kokkos/ModeKokkos_impl.h index 068ab04ce..3a6ae1b81 100644 --- a/src/libnnpif/Kokkos/ModeKokkos_impl.h +++ b/src/libnnpif/Kokkos/ModeKokkos_impl.h @@ -15,6 +15,9 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . +#ifndef MODE_KOKKOS_IMPL_H +#define MODE_KOKKOS_IMPL_H + #include "utility.h" #include // std::min, std::max @@ -31,7 +34,7 @@ namespace nnp // TODO: call base, then copy to Views template -void ModeCabana::setupElementMap() +void ModeKokkos::setupElementMap() { log << "\n"; log << "*** SETUP: ELEMENT MAP ******************" @@ -57,7 +60,7 @@ void ModeCabana::setupElementMap() } template -void ModeCabana::setupElements() +void ModeKokkos::setupElements() { log << "\n"; log << "*** SETUP: ELEMENTS *********************" @@ -75,7 +78,7 @@ void ModeCabana::setupElements() for ( size_t i = 0; i < numElements; ++i ) { - elements.push_back( ElementCabana( i ) ); + elements.push_back( ElementKokkos( i ) ); } if ( settings.keywordExists( "atom_energy" ) ) @@ -110,7 +113,7 @@ void ModeCabana::setupElements() } template -void ModeCabana::setupSymmetryFunctions() +void ModeKokkos::setupSymmetryFunctions() { maxSFperElem = 0; h_numSFperElem = @@ -183,7 +186,7 @@ void ModeCabana::setupSymmetryFunctions() log << "\n"; maxCutoffRadius = 0.0; - for ( vector::iterator it = elements.begin(); it != elements.end(); + for ( vector::iterator it = elements.begin(); it != elements.end(); ++it ) { int attype = it->getIndex(); @@ -234,7 +237,7 @@ void ModeCabana::setupSymmetryFunctions() // TODO: call base, then copy to View template -void ModeCabana::setupSymmetryFunctionScaling( string const &fileName ) +void ModeKokkos::setupSymmetryFunctionScaling( string const &fileName ) { log << "\n"; log << "*** SETUP: SYMMETRY FUNCTION SCALING ****" @@ -341,7 +344,7 @@ void ModeCabana::setupSymmetryFunctionScaling( string const &fileName log << "Smax .... Desired maximum scaled symmetry function value.\n"; log << "t ....... Scaling type.\n"; log << "\n"; - for ( vector::iterator it = elements.begin(); it != elements.end(); + for ( vector::iterator it = elements.begin(); it != elements.end(); ++it ) { int attype = it->getIndex(); @@ -378,7 +381,7 @@ void ModeCabana::setupSymmetryFunctionScaling( string const &fileName } template -void ModeCabana::setupSymmetryFunctionGroups() +void ModeKokkos::setupSymmetryFunctionGroups() { log << "\n"; log << "*** SETUP: SYMMETRY FUNCTION GROUPS *****" @@ -408,7 +411,7 @@ void ModeCabana::setupSymmetryFunctionGroups() h_numSFGperElem = h_t_int( "Mode::numSymmetryFunctionGroupsPerElement", numElements ); - for ( vector::iterator it = elements.begin(); it != elements.end(); + for ( vector::iterator it = elements.begin(); it != elements.end(); ++it ) { int attype = it->getIndex(); @@ -440,7 +443,7 @@ void ModeCabana::setupSymmetryFunctionGroups() } template -void ModeCabana::setupNeuralNetwork() +void ModeKokkos::setupNeuralNetwork() { log << "\n"; log << "*** SETUP: NEURAL NETWORKS **************" @@ -482,7 +485,7 @@ void ModeCabana::setupNeuralNetwork() log << "-----------------------------------------" "--------------------------------------\n"; - for ( vector::iterator it = elements.begin(); it != elements.end(); + for ( vector::iterator it = elements.begin(); it != elements.end(); ++it ) { int attype = it->getIndex(); @@ -527,7 +530,7 @@ void ModeCabana::setupNeuralNetwork() } template -void ModeCabana::setupNeuralNetworkWeights( string const &fileNameFormat ) +void ModeKokkos::setupNeuralNetworkWeights( string const &fileNameFormat ) { log << "\n"; log << "*** SETUP: NEURAL NETWORK WEIGHTS *******" @@ -538,7 +541,7 @@ void ModeCabana::setupNeuralNetworkWeights( string const &fileNameForm fileNameFormat.c_str() ); int count = 0; int AN = 0; - for ( vector::iterator it = elements.begin(); it != elements.end(); + for ( vector::iterator it = elements.begin(); it != elements.end(); ++it ) { const char *estring = elementStrings[count].c_str(); @@ -603,12 +606,12 @@ void ModeCabana::setupNeuralNetworkWeights( string const &fileNameForm template template -void ModeCabana::calculateSymmetryFunctionGroups( + class t_neigh_list> +void ModeKokkos::calculateSymmetryFunctionGroups( t_slice_x x, t_slice_type type, t_slice_G G, t_neigh_list neigh_list, - int N_local, t_neigh_parallel neigh_op_tag, t_angle_parallel angle_op_tag ) + int N_local ) { - Cabana::deep_copy( G, 0.0 ); + Kokkos::deep_copy( G, 0.0 ); Kokkos::RangePolicy policy( 0, N_local ); @@ -668,9 +671,8 @@ void ModeCabana::calculateSymmetryFunctionGroups( } } }; - Cabana::neighbor_parallel_for( - policy, calc_radial_symm_op, neigh_list, Cabana::FirstNeighborsTag(), - neigh_op_tag, "Mode::calculateRadialSymmetryFunctionGroups" ); + Kokkos::parallel_for( + policy, calc_radial_symm_op, "Mode::calculateRadialSymmetryFunctionGroups" ); Kokkos::fence(); auto calc_angular_symm_op = @@ -789,9 +791,8 @@ void ModeCabana::calculateSymmetryFunctionGroups( } } }; - Cabana::neighbor_parallel_for( - policy, calc_angular_symm_op, neigh_list, Cabana::SecondNeighborsTag(), - angle_op_tag, "Mode::calculateAngularSymmetryFunctionGroups" ); + Kokkos::parallel_for( policy, calc_angular_symm_op, + "Mode::calculateAngularSymmetryFunctionGroups" ); Kokkos::fence(); auto scale_symm_op = KOKKOS_LAMBDA( const int i ) @@ -833,7 +834,7 @@ void ModeCabana::calculateSymmetryFunctionGroups( template template -void ModeCabana::calculateAtomicNeuralNetworks( +void ModeKokkos::calculateAtomicNeuralNetworks( t_slice_type type, t_slice_G G, t_slice_dEdG dEdG, t_slice_E E, int N_local ) { auto NN = d_t_NN( "Mode::NN", N_local, numLayers, maxNeurons ); @@ -927,12 +928,10 @@ void ModeCabana::calculateAtomicNeuralNetworks( template template -void ModeCabana::calculateForces( + class t_slice_dEdG, class t_neigh_list> +void ModeKokkos::calculateForces( t_slice_x x, t_slice_f f_a, t_slice_type type, t_slice_dEdG dEdG, - t_neigh_list neigh_list, int N_local, t_neigh_parallel neigh_op_tag, - t_angle_parallel angle_op_tag ) + t_neigh_list neigh_list, int N_local ) { double convForce_ = convLength / convEnergy; @@ -1015,9 +1014,8 @@ void ModeCabana::calculateForces( } } }; - Cabana::neighbor_parallel_for( policy, calc_radial_force_op, neigh_list, - Cabana::FirstNeighborsTag(), neigh_op_tag, - "Mode::calculateRadialForces" ); + Kokkos::parallel_for( policy, calc_radial_force_op, + "Mode::calculateRadialForces" ); Kokkos::fence(); auto calc_angular_force_op = @@ -1208,12 +1206,12 @@ void ModeCabana::calculateForces( } } }; - Cabana::neighbor_parallel_for( policy, calc_angular_force_op, neigh_list, - Cabana::SecondNeighborsTag(), angle_op_tag, - "Mode::calculateAngularForces" ); + Kokkos::parallel_for( policy, calc_angular_force_op, + "Mode::calculateAngularForces" ); Kokkos::fence(); return; } } +#endif diff --git a/src/libnnpif/Kokkos/typesKokkos.h b/src/libnnpif/Kokkos/typesKokkos.h index 66d32bd2d..a177ef023 100644 --- a/src/libnnpif/Kokkos/typesKokkos.h +++ b/src/libnnpif/Kokkos/typesKokkos.h @@ -15,8 +15,8 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . -#ifndef TYPESCABANAMD_H -#define TYPESCABANAMD_H +#ifndef TYPESKOKKOS_H +#define TYPESKOKKOS_H constexpr double CFLENGTH = 1.889726; constexpr double CFENERGY = 0.036749; From 74029cefd238d53b630a6fadbefb915b6eae4c86 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Sun, 31 Jan 2021 00:12:10 -0500 Subject: [PATCH 3/6] Add Kokkos build --- src/libnnpif/makefile | 9 ++++++++- src/makefile | 3 ++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/src/libnnpif/makefile b/src/libnnpif/makefile index 55ec87e6d..3d770eb66 100644 --- a/src/libnnpif/makefile +++ b/src/libnnpif/makefile @@ -31,8 +31,9 @@ include $(PROJECT_DIR)/src/makefile.$(COMP) # makefile or command line. Provide a space-separated list of the following # available options: # * LAMMPS +# * Kokkos # * CabanaMD -INTERFACES=LAMMPS CabanaMD +INTERFACES=LAMMPS Kokkos CabanaMD ############################################################################### @@ -76,6 +77,12 @@ SRC:=$(SRC) $(wildcard $(IF)/*.cpp) HEADERS:=$(HEADERS) $(wildcard $(IF)/*.h) endif +IF=Kokkos +ifneq ($(filter $(IF),$(INTERFACES)),) +SRC:=$(SRC) $(wildcard $(IF)/*.cpp) +HEADERS:=$(HEADERS) $(wildcard $(IF)/*.h) +endif + IF=CabanaMD ifneq ($(filter $(IF),$(INTERFACES)),) SRC:=$(SRC) $(wildcard $(IF)/*.cpp) diff --git a/src/makefile b/src/makefile index 580c91d4e..ff34771cf 100644 --- a/src/makefile +++ b/src/makefile @@ -35,8 +35,9 @@ INSTALL_LIB=$(HOME)/local/lib # Select which interfaces should be compiled into libnnpif. Provide a # space-separated list of the following available options: # * LAMMPS +# * Kokkos # * CabanaMD -INTERFACES=LAMMPS CabanaMD +INTERFACES=LAMMPS Kokkos CabanaMD ############################################################################### From 461b25c71715bf814d5e60adccf3f64a0c4e0a37 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 13 Jan 2021 16:49:39 -0500 Subject: [PATCH 4/6] Add back Cabana specific functions, inheriting from Kokkos version --- src/libnnpif/CabanaMD/ModeCabana.h | 113 +++++ src/libnnpif/CabanaMD/ModeCabana_impl.h | 552 ++++++++++++++++++++++++ 2 files changed, 665 insertions(+) create mode 100644 src/libnnpif/CabanaMD/ModeCabana.h create mode 100644 src/libnnpif/CabanaMD/ModeCabana_impl.h diff --git a/src/libnnpif/CabanaMD/ModeCabana.h b/src/libnnpif/CabanaMD/ModeCabana.h new file mode 100644 index 000000000..23f1e6612 --- /dev/null +++ b/src/libnnpif/CabanaMD/ModeCabana.h @@ -0,0 +1,113 @@ +// n2p2 - A neural network potential package +// Copyright (C) 2018 Andreas Singraber (University of Vienna) +// Copyright (C) 2020 Saaketh Desai and Sam Reeve +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef MODE_CABANA_H +#define MODE_CABANA_H + +#include "ModeKokkos.h" + +#include +#include + +#include "CutoffFunction.h" +#include "Log.h" +#include "Mode.h" +#include "Settings.h" + +#include // std::size_t +#include // std::string +#include // std::vector + +namespace nnp +{ + +/** Derived Cabana main NNP class. + * + * The main n2p2 functions for computing energies and forces are replaced + * to use the Kokkos and Cabana libraries. Only the main compute kernels + * differ from the Kokkos version. + */ +template +class ModeCabana : public ModeKokkos +{ + + public: + using ModeKokkos::ModeKokkos; + using exe_space = typename ModeKokkos::exe_space; + + // Symmetry function Kokkos::Views + using ModeKokkos::d_SF; + using ModeKokkos::d_SFGmemberlist; + using ModeKokkos::d_SFscaling; + + // Per type Kokkos::Views + using ModeKokkos::numSFGperElem; + + using ModeKokkos::cutoffAlpha; + using ModeKokkos::convEnergy; + using ModeKokkos::convLength; + using ModeKokkos::maxSFperElem; + using ModeKokkos::cutoffType; + + /** Calculate forces for all atoms in given structure. + * + * @param[in] x Cabana slice of atom positions. + * @param[in] f Cabana slice of atom forces. + * @param[in] type Cabana slice of atom types. + * @param[in] dEdG Cabana slice of the derivative of energy with respect + to symmetry functions per atom. + * @param[in] N_local Number of atoms. + * @param[in] neigh_op Cabana tag for neighbor parallelism. + * @param[in] angle_op Cabana tag for angular parallelism. + * + * Combine intermediate results from symmetry function and neural network + * computation to atomic forces. Results are stored in f. + */ + template + void calculateForces(t_slice_x x, t_slice_f f, t_slice_type type, + t_slice_dEdG dEdG, t_neigh_list neigh_list, int N_local, + t_neigh_parallel neigh_op, t_angle_parallel angle_op); + + /** Calculate all symmetry function groups for all atoms in given + * structure. + * + * @param[in] x Cabana slice of atom positions. + * @param[in] type Cabana slice of atom types. + * @param[in] G Cabana slice of symmetry functions per atom. + * @param[in] neigh_list neighbor list. + * @param[in] N_local Number of atoms. + * @param[in] neigh_op Cabana tag for neighbor parallelism. + * @param[in] angle_op Cabana tag for angular parallelism. + * + * Note there is no calculateSymmetryFunctions() within this derived class. + * Results are stored in G. + */ + template + void calculateSymmetryFunctionGroups(t_slice_x x, t_slice_type type, + t_slice_G G, t_neigh_list neigh_list, + int N_local, t_neigh_parallel neigh_op, + t_angle_parallel angle_op); +}; + +} + +#include "ModeCabana_impl.h" + +#endif diff --git a/src/libnnpif/CabanaMD/ModeCabana_impl.h b/src/libnnpif/CabanaMD/ModeCabana_impl.h new file mode 100644 index 000000000..87b818af9 --- /dev/null +++ b/src/libnnpif/CabanaMD/ModeCabana_impl.h @@ -0,0 +1,552 @@ +// n2p2 - A neural network potential package +// Copyright (C) 2018 Andreas Singraber (University of Vienna) +// Copyright (C) 2020 Saaketh Desai and Sam Reeve +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#include "utility.h" + +#include // std::min, std::max +#include // atoi, atof +#include // std::ifstream +#include // std::numeric_limits +#include // std::runtime_error +#include + +using namespace std; + +namespace nnp +{ + +template +template +void ModeCabana::calculateSymmetryFunctionGroups( + t_slice_x x, t_slice_type type, t_slice_G G, t_neigh_list neigh_list, + int N_local, t_neigh_parallel neigh_op_tag, t_angle_parallel angle_op_tag ) +{ + Cabana::deep_copy( G, 0.0 ); + + Kokkos::RangePolicy policy( 0, N_local ); + + // Create local copies for lambda + auto numSFGperElem_ = numSFGperElem; + auto SFGmemberlist_ = d_SFGmemberlist; + auto SF_ = d_SF; + auto SFscaling_ = d_SFscaling; + auto maxSFperElem_ = maxSFperElem; + auto convLength_ = convLength; + auto cutoffType_ = cutoffType; + auto cutoffAlpha_ = cutoffAlpha; + + auto calc_radial_symm_op = KOKKOS_LAMBDA( const int i, const int j ) + { + double pfcij = 0.0; + double pdfcij = 0.0; + double eta, rs; + size_t nej; + int memberindex, globalIndex; + double rij, r2ij; + T_F_FLOAT dxij, dyij, dzij; + + int attype = type( i ); + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); + ++groupIndex ) + { + if ( SF_( attype, SFGmemberlist_( attype, groupIndex, 0 ), 1 ) == + 2 ) + { + size_t memberindex0 = SFGmemberlist_( attype, groupIndex, 0 ); + size_t e1 = SF_( attype, memberindex0, 2 ); + double rc = SF_( attype, memberindex0, 7 ); + size_t size = + SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); + + nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( e1 == nej && rij < rc ) + { + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, false ); + for ( size_t k = 0; k < size; ++k ) + { + memberindex = SFGmemberlist_( attype, groupIndex, k ); + globalIndex = SF_( attype, memberindex, 14 ); + eta = SF_( attype, memberindex, 4 ); + rs = SF_( attype, memberindex, 8 ); + G( i, globalIndex ) += + exp( -eta * ( rij - rs ) * ( rij - rs ) ) * pfcij; + } + } + } + } + }; + Cabana::neighbor_parallel_for( + policy, calc_radial_symm_op, neigh_list, Cabana::FirstNeighborsTag(), + neigh_op_tag, "Mode::calculateRadialSymmetryFunctionGroups" ); + Kokkos::fence(); + + auto calc_angular_symm_op = + KOKKOS_LAMBDA( const int i, const int j, const int k ) + { + double pfcij = 0.0, pdfcij = 0.0; + double pfcik = 0.0, pdfcik = 0.0; + double pfcjk = 0.0, pdfcjk = 0.0; + size_t nej, nek; + int memberindex, globalIndex; + double rij, r2ij, rik, r2ik, rjk, r2jk; + T_F_FLOAT dxij, dyij, dzij, dxik, dyik, dzik, dxjk, dyjk, dzjk; + double eta, rs, lambda, zeta; + + int attype = type( i ); + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); + ++groupIndex ) + { + if ( SF_( attype, SFGmemberlist_( attype, groupIndex, 0 ), 1 ) == + 3 ) + { + size_t memberindex0 = SFGmemberlist_( attype, groupIndex, 0 ); + size_t e1 = SF_( attype, memberindex0, 2 ); + size_t e2 = SF_( attype, memberindex0, 3 ); + double rc = SF_( attype, memberindex0, 7 ); + size_t size = + SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); + + nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( ( e1 == nej || e2 == nej ) && rij < rc ) + { + // Calculate cutoff function and derivative. + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, false ); + + nek = type( k ); + + if ( ( e1 == nej && e2 == nek ) || + ( e2 == nej && e1 == nek ) ) + { + dxik = + ( x( i, 0 ) - x( k, 0 ) ) * CFLENGTH * convLength_; + dyik = + ( x( i, 1 ) - x( k, 1 ) ) * CFLENGTH * convLength_; + dzik = + ( x( i, 2 ) - x( k, 2 ) ) * CFLENGTH * convLength_; + r2ik = dxik * dxik + dyik * dyik + dzik * dzik; + rik = sqrt( r2ik ); + if ( rik < rc ) + { + dxjk = dxik - dxij; + dyjk = dyik - dyij; + dzjk = dzik - dzij; + r2jk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; + if ( r2jk < rc * rc ) + { + // Energy calculation. + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcik, pdfcik, rik, rc, false ); + + rjk = sqrt( r2jk ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcjk, pdfcjk, rjk, rc, false ); + + double const rinvijik = 1.0 / rij / rik; + double const costijk = + ( dxij * dxik + dyij * dyik + + dzij * dzik ) * + rinvijik; + double vexp = 0.0, rijs = 0.0, riks = 0.0, + rjks = 0.0; + for ( size_t l = 0; l < size; ++l ) + { + globalIndex = + SF_( attype, + SFGmemberlist_( attype, + groupIndex, l ), + 14 ); + memberindex = SFGmemberlist_( + attype, groupIndex, l ); + eta = SF_( attype, memberindex, 4 ); + lambda = SF_( attype, memberindex, 5 ); + zeta = SF_( attype, memberindex, 6 ); + rs = SF_( attype, memberindex, 8 ); + if ( rs > 0.0 ) + { + rijs = rij - rs; + riks = rik - rs; + rjks = rjk - rs; + vexp = exp( -eta * ( rijs * rijs + + riks * riks + + rjks * rjks ) ); + } + else + vexp = exp( -eta * + ( r2ij + r2ik + r2jk ) ); + double const plambda = + 1.0 + lambda * costijk; + double fg = vexp; + if ( plambda <= 0.0 ) + fg = 0.0; + else + fg *= pow( plambda, ( zeta - 1.0 ) ); + G( i, globalIndex ) += + fg * plambda * pfcij * pfcik * pfcjk; + } // l + } // rjk <= rc + } // rik <= rc + } // elem + } // rij <= rc + } + } + }; + Cabana::neighbor_parallel_for( + policy, calc_angular_symm_op, neigh_list, Cabana::SecondNeighborsTag(), + angle_op_tag, "Mode::calculateAngularSymmetryFunctionGroups" ); + Kokkos::fence(); + + auto scale_symm_op = KOKKOS_LAMBDA( const int i ) + { + int attype = type( i ); + + int memberindex0; + int memberindex, globalIndex; + double raw_value = 0.0; + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); + ++groupIndex ) + { + memberindex0 = SFGmemberlist_( attype, groupIndex, 0 ); + + size_t size = SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); + for ( size_t k = 0; k < size; ++k ) + { + globalIndex = SF_( + attype, SFGmemberlist_( attype, groupIndex, k ), 14 ); + memberindex = SFGmemberlist_( attype, groupIndex, k ); + + if ( SF_( attype, memberindex0, 1 ) == 2 ) + raw_value = G( i, globalIndex ); + else if ( SF_( attype, memberindex0, 1 ) == 3 ) + raw_value = + G( i, globalIndex ) * + pow( 2, ( 1 - SF_( attype, memberindex, 6 ) ) ); + + G( i, globalIndex ) = + this->scale( attype, raw_value, memberindex, SFscaling_ ); + } + } + }; + Kokkos::parallel_for( "Mode::scaleSymmetryFunctionGroups", policy, + scale_symm_op ); + Kokkos::fence(); +} + +template +template +void ModeCabana::calculateForces( + t_slice_x x, t_slice_f f_a, t_slice_type type, t_slice_dEdG dEdG, + t_neigh_list neigh_list, int N_local, t_neigh_parallel neigh_op_tag, + t_angle_parallel angle_op_tag ) +{ + double convForce_ = convLength / convEnergy; + + Kokkos::RangePolicy policy( 0, N_local ); + + // Create local copies for lambda + auto numSFGperElem_ = numSFGperElem; + auto SFGmemberlist_ = d_SFGmemberlist; + auto SF_ = d_SF; + auto SFscaling_ = d_SFscaling; + auto maxSFperElem_ = maxSFperElem; + auto convLength_ = convLength; + auto cutoffType_ = cutoffType; + auto cutoffAlpha_ = cutoffAlpha; + + auto calc_radial_force_op = KOKKOS_LAMBDA( const int i, const int j ) + { + double pfcij = 0.0; + double pdfcij = 0.0; + double rij, r2ij; + T_F_FLOAT dxij, dyij, dzij; + double eta, rs; + int memberindex, globalIndex; + + int attype = type( i ); + + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); + ++groupIndex ) + { + if ( SF_( attype, SFGmemberlist_( attype, groupIndex, 0 ), 1 ) == + 2 ) + { + size_t memberindex0 = SFGmemberlist_( attype, groupIndex, 0 ); + size_t e1 = SF_( attype, memberindex0, 2 ); + double rc = SF_( attype, memberindex0, 7 ); + size_t size = + SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); + + size_t nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( e1 == nej && rij < rc ) + { + // Energy calculation. + // Calculate cutoff function and derivative. + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, true ); + for ( size_t k = 0; k < size; ++k ) + { + globalIndex = SF_( + attype, SFGmemberlist_( attype, groupIndex, k ), + 14 ); + memberindex = SFGmemberlist_( attype, groupIndex, k ); + eta = SF_( attype, memberindex, 4 ); + rs = SF_( attype, memberindex, 8 ); + double pexp = exp( -eta * ( rij - rs ) * ( rij - rs ) ); + // Force calculation. + double const p1 = + SFscaling_( attype, memberindex, 6 ) * + ( pdfcij - 2.0 * eta * ( rij - rs ) * pfcij ) * + pexp / rij; + f_a( i, 0 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dxij ) * CFFORCE * convForce_ ); + f_a( i, 1 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dyij ) * CFFORCE * convForce_ ); + f_a( i, 2 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dzij ) * CFFORCE * convForce_ ); + + f_a( j, 0 ) += ( dEdG( i, globalIndex ) * + ( p1 * dxij ) * CFFORCE * convForce_ ); + f_a( j, 1 ) += ( dEdG( i, globalIndex ) * + ( p1 * dyij ) * CFFORCE * convForce_ ); + f_a( j, 2 ) += ( dEdG( i, globalIndex ) * + ( p1 * dzij ) * CFFORCE * convForce_ ); + } + } + } + } + }; + Cabana::neighbor_parallel_for( policy, calc_radial_force_op, neigh_list, + Cabana::FirstNeighborsTag(), neigh_op_tag, + "Mode::calculateRadialForces" ); + Kokkos::fence(); + + auto calc_angular_force_op = + KOKKOS_LAMBDA( const int i, const int j, const int k ) + { + double pfcij = 0.0; + double pdfcij = 0.0; + double pfcik, pdfcik, pfcjk, pdfcjk; + size_t nej, nek; + double rij, r2ij, rik, r2ik, rjk, r2jk; + T_F_FLOAT dxij, dyij, dzij, dxik, dyik, dzik, dxjk, dyjk, dzjk; + double eta, rs, lambda, zeta; + int memberindex, globalIndex; + + int attype = type( i ); + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); + ++groupIndex ) + { + if ( SF_( attype, SFGmemberlist_( attype, groupIndex, 0 ), 1 ) == + 3 ) + { + size_t memberindex0 = SFGmemberlist_( attype, groupIndex, 0 ); + size_t e1 = SF_( attype, memberindex0, 2 ); + size_t e2 = SF_( attype, memberindex0, 3 ); + double rc = SF_( attype, memberindex0, 7 ); + size_t size = + SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); + + nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( ( e1 == nej || e2 == nej ) && rij < rc ) + { + // Calculate cutoff function and derivative. + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, true ); + + nek = type( k ); + if ( ( e1 == nej && e2 == nek ) || + ( e2 == nej && e1 == nek ) ) + { + dxik = + ( x( i, 0 ) - x( k, 0 ) ) * CFLENGTH * convLength_; + dyik = + ( x( i, 1 ) - x( k, 1 ) ) * CFLENGTH * convLength_; + dzik = + ( x( i, 2 ) - x( k, 2 ) ) * CFLENGTH * convLength_; + r2ik = dxik * dxik + dyik * dyik + dzik * dzik; + rik = sqrt( r2ik ); + if ( rik < rc ) + { + dxjk = dxik - dxij; + dyjk = dyik - dyij; + dzjk = dzik - dzij; + r2jk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; + if ( r2jk < rc * rc ) + { + // Energy calculation. + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcik, pdfcik, rik, rc, true ); + rjk = sqrt( r2jk ); + + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcjk, pdfcjk, rjk, rc, true ); + + double const rinvijik = 1.0 / rij / rik; + double const costijk = + ( dxij * dxik + dyij * dyik + + dzij * dzik ) * + rinvijik; + double const pfc = pfcij * pfcik * pfcjk; + double const r2sum = r2ij + r2ik + r2jk; + double const pr1 = pfcik * pfcjk * pdfcij / rij; + double const pr2 = pfcij * pfcjk * pdfcik / rik; + double const pr3 = pfcij * pfcik * pdfcjk / rjk; + double vexp = 0.0, rijs = 0.0, riks = 0.0, + rjks = 0.0; + for ( size_t l = 0; l < size; ++l ) + { + globalIndex = + SF_( attype, + SFGmemberlist_( attype, + groupIndex, l ), + 14 ); + memberindex = SFGmemberlist_( + attype, groupIndex, l ); + rs = SF_( attype, memberindex, 8 ); + eta = SF_( attype, memberindex, 4 ); + lambda = SF_( attype, memberindex, 5 ); + zeta = SF_( attype, memberindex, 6 ); + if ( rs > 0.0 ) + { + rijs = rij - rs; + riks = rik - rs; + rjks = rjk - rs; + vexp = exp( -eta * ( rijs * rijs + + riks * riks + + rjks * rjks ) ); + } + else + vexp = exp( -eta * r2sum ); + + double const plambda = + 1.0 + lambda * costijk; + double fg = vexp; + if ( plambda <= 0.0 ) + fg = 0.0; + else + fg *= pow( plambda, ( zeta - 1.0 ) ); + + fg *= pow( 2, ( 1 - zeta ) ) * + SFscaling_( attype, memberindex, 6 ); + double const pfczl = pfc * zeta * lambda; + double factorDeriv = + 2.0 * eta / zeta / lambda; + double const p2etapl = + plambda * factorDeriv; + double p1, p2, p3; + if ( rs > 0.0 ) + { + p1 = fg * + ( pfczl * + ( rinvijik - costijk / r2ij - + p2etapl * rijs / rij ) + + pr1 * plambda ); + p2 = fg * + ( pfczl * + ( rinvijik - costijk / r2ik - + p2etapl * riks / rik ) + + pr2 * plambda ); + p3 = + fg * + ( pfczl * ( rinvijik + + p2etapl * rjks / rjk ) - + pr3 * plambda ); + } + else + { + p1 = fg * ( pfczl * ( rinvijik - + costijk / r2ij - + p2etapl ) + + pr1 * plambda ); + p2 = fg * ( pfczl * ( rinvijik - + costijk / r2ik - + p2etapl ) + + pr2 * plambda ); + p3 = fg * + ( pfczl * ( rinvijik + p2etapl ) - + pr3 * plambda ); + } + f_a( i, 0 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dxij + p2 * dxik ) * + CFFORCE * convForce_ ); + f_a( i, 1 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dyij + p2 * dyik ) * + CFFORCE * convForce_ ); + f_a( i, 2 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dzij + p2 * dzik ) * + CFFORCE * convForce_ ); + + f_a( j, 0 ) += ( dEdG( i, globalIndex ) * + ( p1 * dxij + p3 * dxjk ) * + CFFORCE * convForce_ ); + f_a( j, 1 ) += ( dEdG( i, globalIndex ) * + ( p1 * dyij + p3 * dyjk ) * + CFFORCE * convForce_ ); + f_a( j, 2 ) += ( dEdG( i, globalIndex ) * + ( p1 * dzij + p3 * dzjk ) * + CFFORCE * convForce_ ); + + f_a( k, 0 ) += ( dEdG( i, globalIndex ) * + ( p2 * dxik - p3 * dxjk ) * + CFFORCE * convForce_ ); + f_a( k, 1 ) += ( dEdG( i, globalIndex ) * + ( p2 * dyik - p3 * dyjk ) * + CFFORCE * convForce_ ); + f_a( k, 2 ) += ( dEdG( i, globalIndex ) * + ( p2 * dzik - p3 * dzjk ) * + CFFORCE * convForce_ ); + } // l + } // rjk <= rc + } // rik <= rc + } // elem + } // rij <= rc + } + } + }; + Cabana::neighbor_parallel_for( policy, calc_angular_force_op, neigh_list, + Cabana::SecondNeighborsTag(), angle_op_tag, + "Mode::calculateAngularForces" ); + Kokkos::fence(); + + return; +} + +} From b5492ab9ac3436b02f63d754b8a909335688e313 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 19 Feb 2021 20:14:16 -0800 Subject: [PATCH 5/6] Bring up to date with current libnnp --- src/libnnpif/Kokkos/ElementKokkos_impl.h | 1 - src/libnnpif/Kokkos/ModeKokkos.h | 3 ++- src/libnnpif/Kokkos/ModeKokkos_impl.h | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/libnnpif/Kokkos/ElementKokkos_impl.h b/src/libnnpif/Kokkos/ElementKokkos_impl.h index f952b608a..2135bf22f 100644 --- a/src/libnnpif/Kokkos/ElementKokkos_impl.h +++ b/src/libnnpif/Kokkos/ElementKokkos_impl.h @@ -33,7 +33,6 @@ ElementKokkos::ElementKokkos( size_t const _index ) { index = _index; atomicEnergyOffset = 0.0; - neuralNetwork = NULL; } ElementKokkos::~ElementKokkos() {} diff --git a/src/libnnpif/Kokkos/ModeKokkos.h b/src/libnnpif/Kokkos/ModeKokkos.h index 33283e580..0708b084c 100644 --- a/src/libnnpif/Kokkos/ModeKokkos.h +++ b/src/libnnpif/Kokkos/ModeKokkos.h @@ -129,7 +129,8 @@ class ModeKokkos : public Mode * per line, see NeuralNetwork::setConnections() for the correct order. */ void setupNeuralNetworkWeights( - std::string const &fileNameFormat = "weights.%03zu.data") override; + std::string const &fileNameFormatShort = "weights.%03zu.data", + std::string const &fileNameFormatCharge = "weightse.%03zu.data" ) override; /* Compute cutoff for a single atom pair. * diff --git a/src/libnnpif/Kokkos/ModeKokkos_impl.h b/src/libnnpif/Kokkos/ModeKokkos_impl.h index 3a6ae1b81..17e367259 100644 --- a/src/libnnpif/Kokkos/ModeKokkos_impl.h +++ b/src/libnnpif/Kokkos/ModeKokkos_impl.h @@ -530,7 +530,8 @@ void ModeKokkos::setupNeuralNetwork() } template -void ModeKokkos::setupNeuralNetworkWeights( string const &fileNameFormat ) +void ModeKokkos::setupNeuralNetworkWeights( string const &fileNameFormat, + string const &) { log << "\n"; log << "*** SETUP: NEURAL NETWORK WEIGHTS *******" From 0283d5e2427e30a34b5e6440a1f7b54f35168e90 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Wed, 24 Feb 2021 19:32:34 -0800 Subject: [PATCH 6/6] Add back neighbor loops --- src/libnnpif/Kokkos/ElementKokkos_impl.h | 9 +- src/libnnpif/Kokkos/ModeKokkos_impl.h | 657 ++++++++++++----------- 2 files changed, 355 insertions(+), 311 deletions(-) diff --git a/src/libnnpif/Kokkos/ElementKokkos_impl.h b/src/libnnpif/Kokkos/ElementKokkos_impl.h index 2135bf22f..f929e9417 100644 --- a/src/libnnpif/Kokkos/ElementKokkos_impl.h +++ b/src/libnnpif/Kokkos/ElementKokkos_impl.h @@ -15,6 +15,9 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see . +#ifndef ELEMENT_KOKKOS_IMPL_H +#define ELEMENT_KOKKOS_IMPL_H + #include "utility.h" #include // std::sort, std::min, std::max @@ -28,14 +31,14 @@ using namespace std; namespace nnp { -ElementKokkos::ElementKokkos( size_t const _index ) +inline ElementKokkos::ElementKokkos( size_t const _index ) : Element() { index = _index; atomicEnergyOffset = 0.0; } -ElementKokkos::~ElementKokkos() {} +inline ElementKokkos::~ElementKokkos() {} template void ElementKokkos::addSymmetryFunction( string const ¶meters, @@ -490,3 +493,5 @@ void ElementKokkos::updateSymmetryFunctionStatistics */ } + +#endif diff --git a/src/libnnpif/Kokkos/ModeKokkos_impl.h b/src/libnnpif/Kokkos/ModeKokkos_impl.h index 17e367259..fd231d23c 100644 --- a/src/libnnpif/Kokkos/ModeKokkos_impl.h +++ b/src/libnnpif/Kokkos/ModeKokkos_impl.h @@ -606,10 +606,10 @@ void ModeKokkos::setupNeuralNetworkWeights( string const &fileNameForm } template -template void ModeKokkos::calculateSymmetryFunctionGroups( - t_slice_x x, t_slice_type type, t_slice_G G, t_neigh_list neigh_list, + t_x x, t_type type, t_G G, t_neigh_list neigh_list, int N_local ) { Kokkos::deep_copy( G, 0.0 ); @@ -626,7 +626,7 @@ void ModeKokkos::calculateSymmetryFunctionGroups( auto cutoffType_ = cutoffType; auto cutoffAlpha_ = cutoffAlpha; - auto calc_radial_symm_op = KOKKOS_LAMBDA( const int i, const int j ) + auto calc_radial_symm_op = KOKKOS_LAMBDA( const int i ) { double pfcij = 0.0; double pdfcij = 0.0; @@ -636,7 +636,10 @@ void ModeKokkos::calculateSymmetryFunctionGroups( double rij, r2ij; T_F_FLOAT dxij, dyij, dzij; + typename t_neigh_list::t_neighs neighs_i = neigh_list.get_neighs(i); + const int num_neighs_i = neighs_i.get_num_neighs(); int attype = type( i ); + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); ++groupIndex ) { @@ -649,35 +652,38 @@ void ModeKokkos::calculateSymmetryFunctionGroups( size_t size = SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); - nej = type( j ); - dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; - dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; - dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; - r2ij = dxij * dxij + dyij * dyij + dzij * dzij; - rij = sqrt( r2ij ); - if ( e1 == nej && rij < rc ) + for(int jj = 0; jj < num_neighs_i; jj++) { - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, false ); - for ( size_t k = 0; k < size; ++k ) + T_INT j = neighs_i(jj); + nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( e1 == nej && rij < rc ) { - memberindex = SFGmemberlist_( attype, groupIndex, k ); - globalIndex = SF_( attype, memberindex, 14 ); - eta = SF_( attype, memberindex, 4 ); - rs = SF_( attype, memberindex, 8 ); - G( i, globalIndex ) += - exp( -eta * ( rij - rs ) * ( rij - rs ) ) * pfcij; + compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, false ); + for ( size_t k = 0; k < size; ++k ) + { + memberindex = SFGmemberlist_( attype, groupIndex, k ); + globalIndex = SF_( attype, memberindex, 14 ); + eta = SF_( attype, memberindex, 4 ); + rs = SF_( attype, memberindex, 8 ); + G( i, globalIndex ) += + exp( -eta * ( rij - rs ) * ( rij - rs ) ) * pfcij; + } } } } } }; - Kokkos::parallel_for( - policy, calc_radial_symm_op, "Mode::calculateRadialSymmetryFunctionGroups" ); + Kokkos::parallel_for( "Mode::calculateRadialSymmetryFunctionGroups", + policy, calc_radial_symm_op ); Kokkos::fence(); - auto calc_angular_symm_op = - KOKKOS_LAMBDA( const int i, const int j, const int k ) + auto calc_angular_symm_op = KOKKOS_LAMBDA( const int i ) { double pfcij = 0.0, pdfcij = 0.0; double pfcik = 0.0, pdfcik = 0.0; @@ -688,7 +694,10 @@ void ModeKokkos::calculateSymmetryFunctionGroups( T_F_FLOAT dxij, dyij, dzij, dxik, dyik, dzik, dxjk, dyjk, dzjk; double eta, rs, lambda, zeta; + typename t_neigh_list::t_neighs neighs_i = neigh_list.get_neighs(i); + const int num_neighs_i = neighs_i.get_num_neighs(); int attype = type( i ); + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); ++groupIndex ) { @@ -702,98 +711,108 @@ void ModeKokkos::calculateSymmetryFunctionGroups( size_t size = SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); - nej = type( j ); - dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; - dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; - dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; - r2ij = dxij * dxij + dyij * dyij + dzij * dzij; - rij = sqrt( r2ij ); - if ( ( e1 == nej || e2 == nej ) && rij < rc ) + for(int jj = 0; jj < num_neighs_i; jj++) { - // Calculate cutoff function and derivative. - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, false ); - - nek = type( k ); - - if ( ( e1 == nej && e2 == nek ) || - ( e2 == nej && e1 == nek ) ) + T_INT j = neighs_i(jj); + nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( ( e1 == nej || e2 == nej ) && rij < rc ) { - dxik = - ( x( i, 0 ) - x( k, 0 ) ) * CFLENGTH * convLength_; - dyik = - ( x( i, 1 ) - x( k, 1 ) ) * CFLENGTH * convLength_; - dzik = - ( x( i, 2 ) - x( k, 2 ) ) * CFLENGTH * convLength_; - r2ik = dxik * dxik + dyik * dyik + dzik * dzik; - rik = sqrt( r2ik ); - if ( rik < rc ) + // Calculate cutoff function and derivative. + compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, false ); + + typename t_neigh_list::t_neighs neighs_j = neigh_list.get_neighs(j); + const int num_neighs_j = neighs_j.get_num_neighs(); + for(int kk = 0; kk < num_neighs_j; kk++) { - dxjk = dxik - dxij; - dyjk = dyik - dyij; - dzjk = dzik - dzij; - r2jk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; - if ( r2jk < rc * rc ) + T_INT k = neighs_i(kk); + nek = type( k ); + + if ( ( e1 == nej && e2 == nek ) || + ( e2 == nej && e1 == nek ) ) { - // Energy calculation. - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcik, pdfcik, rik, rc, false ); - - rjk = sqrt( r2jk ); - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcjk, pdfcjk, rjk, rc, false ); - - double const rinvijik = 1.0 / rij / rik; - double const costijk = - ( dxij * dxik + dyij * dyik + - dzij * dzik ) * - rinvijik; - double vexp = 0.0, rijs = 0.0, riks = 0.0, - rjks = 0.0; - for ( size_t l = 0; l < size; ++l ) + dxik = + ( x( i, 0 ) - x( k, 0 ) ) * CFLENGTH * convLength_; + dyik = + ( x( i, 1 ) - x( k, 1 ) ) * CFLENGTH * convLength_; + dzik = + ( x( i, 2 ) - x( k, 2 ) ) * CFLENGTH * convLength_; + r2ik = dxik * dxik + dyik * dyik + dzik * dzik; + rik = sqrt( r2ik ); + if ( rik < rc ) { - globalIndex = - SF_( attype, - SFGmemberlist_( attype, - groupIndex, l ), - 14 ); - memberindex = SFGmemberlist_( - attype, groupIndex, l ); - eta = SF_( attype, memberindex, 4 ); - lambda = SF_( attype, memberindex, 5 ); - zeta = SF_( attype, memberindex, 6 ); - rs = SF_( attype, memberindex, 8 ); - if ( rs > 0.0 ) + dxjk = dxik - dxij; + dyjk = dyik - dyij; + dzjk = dzik - dzij; + r2jk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; + if ( r2jk < rc * rc ) { - rijs = rij - rs; - riks = rik - rs; - rjks = rjk - rs; - vexp = exp( -eta * ( rijs * rijs + - riks * riks + - rjks * rjks ) ); - } - else - vexp = exp( -eta * - ( r2ij + r2ik + r2jk ) ); - double const plambda = - 1.0 + lambda * costijk; - double fg = vexp; - if ( plambda <= 0.0 ) - fg = 0.0; - else - fg *= pow( plambda, ( zeta - 1.0 ) ); - G( i, globalIndex ) += - fg * plambda * pfcij * pfcik * pfcjk; - } // l - } // rjk <= rc - } // rik <= rc - } // elem - } // rij <= rc + // Energy calculation. + compute_cutoff( cutoffType_, cutoffAlpha_, + pfcik, pdfcik, rik, rc, false ); + + rjk = sqrt( r2jk ); + compute_cutoff( cutoffType_, cutoffAlpha_, + pfcjk, pdfcjk, rjk, rc, false ); + + double const rinvijik = 1.0 / rij / rik; + double const costijk = + ( dxij * dxik + dyij * dyik + + dzij * dzik ) * + rinvijik; + double vexp = 0.0, rijs = 0.0, riks = 0.0, + rjks = 0.0; + for ( size_t l = 0; l < size; ++l ) + { + globalIndex = + SF_( attype, + SFGmemberlist_( attype, + groupIndex, l ), + 14 ); + memberindex = SFGmemberlist_( + attype, groupIndex, l ); + eta = SF_( attype, memberindex, 4 ); + lambda = SF_( attype, memberindex, 5 ); + zeta = SF_( attype, memberindex, 6 ); + rs = SF_( attype, memberindex, 8 ); + if ( rs > 0.0 ) + { + rijs = rij - rs; + riks = rik - rs; + rjks = rjk - rs; + vexp = exp( -eta * ( rijs * rijs + + riks * riks + + rjks * rjks ) ); + } + else + vexp = exp( -eta * + ( r2ij + r2ik + r2jk ) ); + double const plambda = + 1.0 + lambda * costijk; + double fg = vexp; + if ( plambda <= 0.0 ) + fg = 0.0; + else + fg *= pow( plambda, ( zeta - 1.0 ) ); + G( i, globalIndex ) += + fg * plambda * pfcij * pfcik * pfcjk; + } // l + } // rjk <= rc + } // rik <= rc + } // elem + } // k + } // rij <= rc + } // j } } }; - Kokkos::parallel_for( policy, calc_angular_symm_op, - "Mode::calculateAngularSymmetryFunctionGroups" ); + Kokkos::parallel_for( "Mode::calculateAngularSymmetryFunctionGroups", + policy, calc_angular_symm_op ); Kokkos::fence(); auto scale_symm_op = KOKKOS_LAMBDA( const int i ) @@ -833,10 +852,10 @@ void ModeKokkos::calculateSymmetryFunctionGroups( } template -template +template void ModeKokkos::calculateAtomicNeuralNetworks( - t_slice_type type, t_slice_G G, t_slice_dEdG dEdG, t_slice_E E, int N_local ) + t_type type, t_G G, t_dEdG dEdG, t_E E, int N_local ) { auto NN = d_t_NN( "Mode::NN", N_local, numLayers, maxNeurons ); auto dfdx = d_t_NN( "Mode::dfdx", N_local, numLayers, maxNeurons ); @@ -928,10 +947,10 @@ void ModeKokkos::calculateAtomicNeuralNetworks( } template -template +template void ModeKokkos::calculateForces( - t_slice_x x, t_slice_f f_a, t_slice_type type, t_slice_dEdG dEdG, + t_x x, t_f f_a, t_type type, t_dEdG dEdG, t_neigh_list neigh_list, int N_local ) { double convForce_ = convLength / convEnergy; @@ -948,7 +967,7 @@ void ModeKokkos::calculateForces( auto cutoffType_ = cutoffType; auto cutoffAlpha_ = cutoffAlpha; - auto calc_radial_force_op = KOKKOS_LAMBDA( const int i, const int j ) + auto calc_radial_force_op = KOKKOS_LAMBDA( const int i ) { double pfcij = 0.0; double pdfcij = 0.0; @@ -957,6 +976,9 @@ void ModeKokkos::calculateForces( double eta, rs; int memberindex, globalIndex; + + typename t_neigh_list::t_neighs neighs_i = neigh_list.get_neighs(i); + const int num_neighs_i = neighs_i.get_num_neighs(); int attype = type( i ); for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); @@ -971,56 +993,59 @@ void ModeKokkos::calculateForces( size_t size = SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); - size_t nej = type( j ); - dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; - dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; - dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; - r2ij = dxij * dxij + dyij * dyij + dzij * dzij; - rij = sqrt( r2ij ); - if ( e1 == nej && rij < rc ) + for(int jj = 0; jj < num_neighs_i; jj++) { - // Energy calculation. - // Calculate cutoff function and derivative. - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, true ); - for ( size_t k = 0; k < size; ++k ) + T_INT j = neighs_i(jj); + size_t nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( e1 == nej && rij < rc ) { - globalIndex = SF_( - attype, SFGmemberlist_( attype, groupIndex, k ), - 14 ); - memberindex = SFGmemberlist_( attype, groupIndex, k ); - eta = SF_( attype, memberindex, 4 ); - rs = SF_( attype, memberindex, 8 ); - double pexp = exp( -eta * ( rij - rs ) * ( rij - rs ) ); - // Force calculation. - double const p1 = - SFscaling_( attype, memberindex, 6 ) * - ( pdfcij - 2.0 * eta * ( rij - rs ) * pfcij ) * - pexp / rij; - f_a( i, 0 ) -= ( dEdG( i, globalIndex ) * - ( p1 * dxij ) * CFFORCE * convForce_ ); - f_a( i, 1 ) -= ( dEdG( i, globalIndex ) * - ( p1 * dyij ) * CFFORCE * convForce_ ); - f_a( i, 2 ) -= ( dEdG( i, globalIndex ) * - ( p1 * dzij ) * CFFORCE * convForce_ ); - - f_a( j, 0 ) += ( dEdG( i, globalIndex ) * - ( p1 * dxij ) * CFFORCE * convForce_ ); - f_a( j, 1 ) += ( dEdG( i, globalIndex ) * - ( p1 * dyij ) * CFFORCE * convForce_ ); - f_a( j, 2 ) += ( dEdG( i, globalIndex ) * - ( p1 * dzij ) * CFFORCE * convForce_ ); + // Energy calculation. + // Calculate cutoff function and derivative. + compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, true ); + for ( size_t k = 0; k < size; ++k ) + { + globalIndex = SF_( + attype, SFGmemberlist_( attype, groupIndex, k ), + 14 ); + memberindex = SFGmemberlist_( attype, groupIndex, k ); + eta = SF_( attype, memberindex, 4 ); + rs = SF_( attype, memberindex, 8 ); + double pexp = exp( -eta * ( rij - rs ) * ( rij - rs ) ); + // Force calculation. + double const p1 = + SFscaling_( attype, memberindex, 6 ) * + ( pdfcij - 2.0 * eta * ( rij - rs ) * pfcij ) * + pexp / rij; + f_a( i, 0 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dxij ) * CFFORCE * convForce_ ); + f_a( i, 1 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dyij ) * CFFORCE * convForce_ ); + f_a( i, 2 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dzij ) * CFFORCE * convForce_ ); + + f_a( j, 0 ) += ( dEdG( i, globalIndex ) * + ( p1 * dxij ) * CFFORCE * convForce_ ); + f_a( j, 1 ) += ( dEdG( i, globalIndex ) * + ( p1 * dyij ) * CFFORCE * convForce_ ); + f_a( j, 2 ) += ( dEdG( i, globalIndex ) * + ( p1 * dzij ) * CFFORCE * convForce_ ); + } } } } } }; - Kokkos::parallel_for( policy, calc_radial_force_op, - "Mode::calculateRadialForces" ); + Kokkos::parallel_for( "Mode::calculateRadialForces", + policy, calc_radial_force_op ); Kokkos::fence(); - auto calc_angular_force_op = - KOKKOS_LAMBDA( const int i, const int j, const int k ) + auto calc_angular_force_op = KOKKOS_LAMBDA( const int i ) { double pfcij = 0.0; double pdfcij = 0.0; @@ -1031,7 +1056,10 @@ void ModeKokkos::calculateForces( double eta, rs, lambda, zeta; int memberindex, globalIndex; + typename t_neigh_list::t_neighs neighs_i = neigh_list.get_neighs(i); + const int num_neighs_i = neighs_i.get_num_neighs(); int attype = type( i ); + for ( int groupIndex = 0; groupIndex < numSFGperElem_( attype ); ++groupIndex ) { @@ -1045,170 +1073,181 @@ void ModeKokkos::calculateForces( size_t size = SFGmemberlist_( attype, groupIndex, maxSFperElem_ ); - nej = type( j ); - dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; - dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; - dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; - r2ij = dxij * dxij + dyij * dyij + dzij * dzij; - rij = sqrt( r2ij ); - if ( ( e1 == nej || e2 == nej ) && rij < rc ) + for(int jj = 0; jj < num_neighs_i; jj++) { - // Calculate cutoff function and derivative. - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, true ); - - nek = type( k ); - if ( ( e1 == nej && e2 == nek ) || - ( e2 == nej && e1 == nek ) ) + T_INT j = neighs_i(jj); + nej = type( j ); + dxij = ( x( i, 0 ) - x( j, 0 ) ) * CFLENGTH * convLength_; + dyij = ( x( i, 1 ) - x( j, 1 ) ) * CFLENGTH * convLength_; + dzij = ( x( i, 2 ) - x( j, 2 ) ) * CFLENGTH * convLength_; + r2ij = dxij * dxij + dyij * dyij + dzij * dzij; + rij = sqrt( r2ij ); + if ( ( e1 == nej || e2 == nej ) && rij < rc ) { - dxik = - ( x( i, 0 ) - x( k, 0 ) ) * CFLENGTH * convLength_; - dyik = - ( x( i, 1 ) - x( k, 1 ) ) * CFLENGTH * convLength_; - dzik = - ( x( i, 2 ) - x( k, 2 ) ) * CFLENGTH * convLength_; - r2ik = dxik * dxik + dyik * dyik + dzik * dzik; - rik = sqrt( r2ik ); - if ( rik < rc ) + // Calculate cutoff function and derivative. + compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, true ); + + typename t_neigh_list::t_neighs neighs_j = neigh_list.get_neighs(j); + const int num_neighs_j = neighs_j.get_num_neighs(); + + for(int kk = 0; kk < num_neighs_j; kk++) { - dxjk = dxik - dxij; - dyjk = dyik - dyij; - dzjk = dzik - dzij; - r2jk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; - if ( r2jk < rc * rc ) + T_INT k = neighs_j(kk); + nek = type( k ); + if ( ( e1 == nej && e2 == nek ) || + ( e2 == nej && e1 == nek ) ) { - // Energy calculation. - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcik, pdfcik, rik, rc, true ); - rjk = sqrt( r2jk ); - - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcjk, pdfcjk, rjk, rc, true ); - - double const rinvijik = 1.0 / rij / rik; - double const costijk = - ( dxij * dxik + dyij * dyik + - dzij * dzik ) * - rinvijik; - double const pfc = pfcij * pfcik * pfcjk; - double const r2sum = r2ij + r2ik + r2jk; - double const pr1 = pfcik * pfcjk * pdfcij / rij; - double const pr2 = pfcij * pfcjk * pdfcik / rik; - double const pr3 = pfcij * pfcik * pdfcjk / rjk; - double vexp = 0.0, rijs = 0.0, riks = 0.0, - rjks = 0.0; - for ( size_t l = 0; l < size; ++l ) + dxik = + ( x( i, 0 ) - x( k, 0 ) ) * CFLENGTH * convLength_; + dyik = + ( x( i, 1 ) - x( k, 1 ) ) * CFLENGTH * convLength_; + dzik = + ( x( i, 2 ) - x( k, 2 ) ) * CFLENGTH * convLength_; + r2ik = dxik * dxik + dyik * dyik + dzik * dzik; + rik = sqrt( r2ik ); + if ( rik < rc ) { - globalIndex = - SF_( attype, - SFGmemberlist_( attype, - groupIndex, l ), - 14 ); - memberindex = SFGmemberlist_( - attype, groupIndex, l ); - rs = SF_( attype, memberindex, 8 ); - eta = SF_( attype, memberindex, 4 ); - lambda = SF_( attype, memberindex, 5 ); - zeta = SF_( attype, memberindex, 6 ); - if ( rs > 0.0 ) - { - rijs = rij - rs; - riks = rik - rs; - rjks = rjk - rs; - vexp = exp( -eta * ( rijs * rijs + - riks * riks + - rjks * rjks ) ); - } - else - vexp = exp( -eta * r2sum ); - - double const plambda = - 1.0 + lambda * costijk; - double fg = vexp; - if ( plambda <= 0.0 ) - fg = 0.0; - else - fg *= pow( plambda, ( zeta - 1.0 ) ); - - fg *= pow( 2, ( 1 - zeta ) ) * - SFscaling_( attype, memberindex, 6 ); - double const pfczl = pfc * zeta * lambda; - double factorDeriv = - 2.0 * eta / zeta / lambda; - double const p2etapl = - plambda * factorDeriv; - double p1, p2, p3; - if ( rs > 0.0 ) - { - p1 = fg * - ( pfczl * - ( rinvijik - costijk / r2ij - - p2etapl * rijs / rij ) + - pr1 * plambda ); - p2 = fg * - ( pfczl * - ( rinvijik - costijk / r2ik - - p2etapl * riks / rik ) + - pr2 * plambda ); - p3 = - fg * - ( pfczl * ( rinvijik + - p2etapl * rjks / rjk ) - - pr3 * plambda ); - } - else + dxjk = dxik - dxij; + dyjk = dyik - dyij; + dzjk = dzik - dzij; + r2jk = dxjk * dxjk + dyjk * dyjk + dzjk * dzjk; + if ( r2jk < rc * rc ) { - p1 = fg * ( pfczl * ( rinvijik - - costijk / r2ij - - p2etapl ) + - pr1 * plambda ); - p2 = fg * ( pfczl * ( rinvijik - - costijk / r2ik - - p2etapl ) + - pr2 * plambda ); - p3 = fg * - ( pfczl * ( rinvijik + p2etapl ) - - pr3 * plambda ); - } - f_a( i, 0 ) -= ( dEdG( i, globalIndex ) * - ( p1 * dxij + p2 * dxik ) * - CFFORCE * convForce_ ); - f_a( i, 1 ) -= ( dEdG( i, globalIndex ) * - ( p1 * dyij + p2 * dyik ) * - CFFORCE * convForce_ ); - f_a( i, 2 ) -= ( dEdG( i, globalIndex ) * - ( p1 * dzij + p2 * dzik ) * - CFFORCE * convForce_ ); - - f_a( j, 0 ) += ( dEdG( i, globalIndex ) * - ( p1 * dxij + p3 * dxjk ) * - CFFORCE * convForce_ ); - f_a( j, 1 ) += ( dEdG( i, globalIndex ) * - ( p1 * dyij + p3 * dyjk ) * - CFFORCE * convForce_ ); - f_a( j, 2 ) += ( dEdG( i, globalIndex ) * - ( p1 * dzij + p3 * dzjk ) * - CFFORCE * convForce_ ); - - f_a( k, 0 ) += ( dEdG( i, globalIndex ) * - ( p2 * dxik - p3 * dxjk ) * - CFFORCE * convForce_ ); - f_a( k, 1 ) += ( dEdG( i, globalIndex ) * - ( p2 * dyik - p3 * dyjk ) * - CFFORCE * convForce_ ); - f_a( k, 2 ) += ( dEdG( i, globalIndex ) * - ( p2 * dzik - p3 * dzjk ) * - CFFORCE * convForce_ ); - } // l - } // rjk <= rc - } // rik <= rc - } // elem - } // rij <= rc + // Energy calculation. + compute_cutoff( cutoffType_, cutoffAlpha_, + pfcik, pdfcik, rik, rc, true ); + rjk = sqrt( r2jk ); + + compute_cutoff( cutoffType_, cutoffAlpha_, + pfcjk, pdfcjk, rjk, rc, true ); + + double const rinvijik = 1.0 / rij / rik; + double const costijk = + ( dxij * dxik + dyij * dyik + + dzij * dzik ) * + rinvijik; + double const pfc = pfcij * pfcik * pfcjk; + double const r2sum = r2ij + r2ik + r2jk; + double const pr1 = pfcik * pfcjk * pdfcij / rij; + double const pr2 = pfcij * pfcjk * pdfcik / rik; + double const pr3 = pfcij * pfcik * pdfcjk / rjk; + double vexp = 0.0, rijs = 0.0, riks = 0.0, + rjks = 0.0; + for ( size_t l = 0; l < size; ++l ) + { + globalIndex = + SF_( attype, + SFGmemberlist_( attype, + groupIndex, l ), + 14 ); + memberindex = SFGmemberlist_( + attype, groupIndex, l ); + rs = SF_( attype, memberindex, 8 ); + eta = SF_( attype, memberindex, 4 ); + lambda = SF_( attype, memberindex, 5 ); + zeta = SF_( attype, memberindex, 6 ); + if ( rs > 0.0 ) + { + rijs = rij - rs; + riks = rik - rs; + rjks = rjk - rs; + vexp = exp( -eta * ( rijs * rijs + + riks * riks + + rjks * rjks ) ); + } + else + vexp = exp( -eta * r2sum ); + + double const plambda = + 1.0 + lambda * costijk; + double fg = vexp; + if ( plambda <= 0.0 ) + fg = 0.0; + else + fg *= pow( plambda, ( zeta - 1.0 ) ); + + fg *= pow( 2, ( 1 - zeta ) ) * + SFscaling_( attype, memberindex, 6 ); + double const pfczl = pfc * zeta * lambda; + double factorDeriv = + 2.0 * eta / zeta / lambda; + double const p2etapl = + plambda * factorDeriv; + double p1, p2, p3; + if ( rs > 0.0 ) + { + p1 = fg * + ( pfczl * + ( rinvijik - costijk / r2ij - + p2etapl * rijs / rij ) + + pr1 * plambda ); + p2 = fg * + ( pfczl * + ( rinvijik - costijk / r2ik - + p2etapl * riks / rik ) + + pr2 * plambda ); + p3 = + fg * + ( pfczl * ( rinvijik + + p2etapl * rjks / rjk ) - + pr3 * plambda ); + } + else + { + p1 = fg * ( pfczl * ( rinvijik - + costijk / r2ij - + p2etapl ) + + pr1 * plambda ); + p2 = fg * ( pfczl * ( rinvijik - + costijk / r2ik - + p2etapl ) + + pr2 * plambda ); + p3 = fg * + ( pfczl * ( rinvijik + p2etapl ) - + pr3 * plambda ); + } + f_a( i, 0 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dxij + p2 * dxik ) * + CFFORCE * convForce_ ); + f_a( i, 1 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dyij + p2 * dyik ) * + CFFORCE * convForce_ ); + f_a( i, 2 ) -= ( dEdG( i, globalIndex ) * + ( p1 * dzij + p2 * dzik ) * + CFFORCE * convForce_ ); + + f_a( j, 0 ) += ( dEdG( i, globalIndex ) * + ( p1 * dxij + p3 * dxjk ) * + CFFORCE * convForce_ ); + f_a( j, 1 ) += ( dEdG( i, globalIndex ) * + ( p1 * dyij + p3 * dyjk ) * + CFFORCE * convForce_ ); + f_a( j, 2 ) += ( dEdG( i, globalIndex ) * + ( p1 * dzij + p3 * dzjk ) * + CFFORCE * convForce_ ); + + f_a( k, 0 ) += ( dEdG( i, globalIndex ) * + ( p2 * dxik - p3 * dxjk ) * + CFFORCE * convForce_ ); + f_a( k, 1 ) += ( dEdG( i, globalIndex ) * + ( p2 * dyik - p3 * dyjk ) * + CFFORCE * convForce_ ); + f_a( k, 2 ) += ( dEdG( i, globalIndex ) * + ( p2 * dzik - p3 * dzjk ) * + CFFORCE * convForce_ ); + } // l + } // rjk <= rc + } // rik <= rc + } // elem + } // k + } // rij <= rc + } // j } } }; - Kokkos::parallel_for( policy, calc_angular_force_op, - "Mode::calculateAngularForces" ); + Kokkos::parallel_for( "Mode::calculateAngularForces", + policy, calc_angular_force_op ); Kokkos::fence(); return;