diff --git a/src/libnnpif/CabanaMD/ModeCabana.h b/src/libnnpif/CabanaMD/ModeCabana.h index 2374c8e01..23f1e6612 100644 --- a/src/libnnpif/CabanaMD/ModeCabana.h +++ b/src/libnnpif/CabanaMD/ModeCabana.h @@ -18,11 +18,10 @@ #ifndef MODE_CABANA_H #define MODE_CABANA_H -#include "./ElementCabana.h" -#include "typesCabana.h" +#include "ModeKokkos.h" -#include #include +#include #include "CutoffFunction.h" #include "Log.h" @@ -39,126 +38,30 @@ 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. Most setup functions are - * overridden; some are replaced as needed to work within device kernels. + * to use the Kokkos and Cabana libraries. Only the main compute kernels + * differ from the Kokkos version. */ template -class ModeCabana : public Mode +class ModeCabana : public ModeKokkos { public: - using Mode::Mode; - - // Kokkos settings - using device_type = t_device; - using exe_space = typename device_type::execution_space; - using memory_space = typename device_type::memory_space; - typedef typename exe_space::array_layout layout; - using host_space = Kokkos::HostSpace; - - // Per-type Kokkos::Views - using d_t_mass = Kokkos::View; - using h_t_mass = Kokkos::View; - using d_t_int = Kokkos::View; - using h_t_int = Kokkos::View; + using ModeKokkos::ModeKokkos; + using exe_space = typename ModeKokkos::exe_space; - // SymmetryFunctionTypes Kokkos::Views - using d_t_SF = Kokkos::View; - using t_SF = Kokkos::View; - using d_t_SFscaling = Kokkos::View; - using t_SFscaling = Kokkos::View; - using d_t_SFGmemberlist = Kokkos::View; - using t_SFGmemberlist = Kokkos::View; - - // NN Kokkos::Views - using d_t_bias = Kokkos::View; - using t_bias = Kokkos::View; - using d_t_weights = Kokkos::View; - using t_weights = Kokkos::View; - using d_t_NN = Kokkos::View; + // Symmetry function Kokkos::Views + using ModeKokkos::d_SF; + using ModeKokkos::d_SFGmemberlist; + using ModeKokkos::d_SFscaling; - /** Set up the element map. - * - * Uses keyword `elements`. This function should follow immediately after - * settings are loaded via loadSettingsFile(). - */ - void setupElementMap() override; - /** Set up all Element instances. - * - * Uses keywords `number_of_elements` and `atom_energy`. This function - * should follow immediately after setupElementMap(). - */ - void setupElements() override; - /** Set up all symmetry functions. - * - * Uses keyword `symfunction_short`. Reads all symmetry functions from - * settings and automatically assigns them to the correct element. - */ - void setupSymmetryFunctions() override; - /** Set up symmetry function scaling from file. - * - * @param[in] fileName Scaling file name. - * - * Uses keywords `scale_symmetry_functions`, `center_symmetry_functions`, - * `scale_symmetry_functions_sigma`, `scale_min_short` and - * `scale_max_short`. Reads in scaling information and sets correct scaling - * behavior for all symmetry functions. Call after - * setupSymmetryFunctions(). - */ - void setupSymmetryFunctionScaling( - std::string const &fileName = "scaling.data") override; - /** Set up symmetry function groups. - * - * Does not use any keywords. Call after setupSymmetryFunctions() and - * ensure that correct scaling behavior has already been set. - */ - void setupSymmetryFunctionGroups() override; - /** Set up neural networks for all elements. - * - * Uses keywords `global_hidden_layers_short`, `global_nodes_short`, - * `global_activation_short`, `normalize_nodes`. Call after - * setupSymmetryFunctions(), only then the number of input layer neurons is - * known. - */ - void setupNeuralNetwork() override; - /** Set up neural network weights from files. - * - * @param[in] fileNameFormat Format for weights file name. The string must - * contain one placeholder for the atomic number. - * - * Does not use any keywords. The weight files should contain one weight - * per line, see NeuralNetwork::setConnections() for the correct order. - */ - void setupNeuralNetworkWeights( - std::string const &fileNameFormat = "weights.%03zu.data") override; + // Per type Kokkos::Views + using ModeKokkos::numSFGperElem; - /* Compute cutoff for a single atom pair. - * - * @param[in] cutoffType Cutoff type. - * @param[in] fc Cutoff function value. - * @param[in] dfc Derivative of cutoff function. - * @param[in] r Atom separation distance. - * @param[in] rc Cutoff radius. - * @param[in] derivative If the cutoff derivative should be computed. - * - * Callable within device kernel. - */ - KOKKOS_INLINE_FUNCTION - void compute_cutoff(CutoffFunction::CutoffType cutoffType, double cutoffAlpha, - double &fc, double &dfc, double r, double rc, - bool derivative) const; - /* - * Scale one symmetry function. - * - * @param[in] attype Atom type. - * @param[in] value Scaled symmetry function value. - * @param[in] k Symmetry function index. - * @param[in] SFscaling Kokkos host View of symmetry function scaling. - * - * Callable within device kernel. - */ - KOKKOS_INLINE_FUNCTION - double scale(int attype, double value, int k, d_t_SFscaling SFscaling) const; + using ModeKokkos::cutoffAlpha; + using ModeKokkos::convEnergy; + using ModeKokkos::convLength; + using ModeKokkos::maxSFperElem; + using ModeKokkos::cutoffType; /** Calculate forces for all atoms in given structure. * @@ -181,30 +84,13 @@ class ModeCabana : public Mode t_slice_dEdG dEdG, t_neigh_list neigh_list, int N_local, t_neigh_parallel neigh_op, t_angle_parallel angle_op); - /** 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 - to symmetry functions per atom. - * @param[in] E Cabana slice of energy per atom. - * @param[in] N_local Number of atoms. - * - * The atomic energy is stored in E. - */ - template - void calculateAtomicNeuralNetworks(t_slice_type type, t_slice_G G, - t_slice_dEdG dEdG, t_slice_E E, - int N_local); - /** 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] 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. @@ -218,128 +104,8 @@ class ModeCabana : public Mode t_slice_G G, t_neigh_list neigh_list, int N_local, t_neigh_parallel neigh_op, t_angle_parallel angle_op); - - using Mode::log; - - /// list of element symbols in order of periodic table - // (duplicated from ElementMap) - std::vector knownElements = { - "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", - "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", - "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", - "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", - "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", - "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", - "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", - "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", - "Bk", "Cf", "Es", "Fm", "Md", "No"}; - - // Symmetry function Kokkos::Views - d_t_SF d_SF; - t_SF SF; - d_t_SFGmemberlist d_SFGmemberlist; - t_SFGmemberlist SFGmemberlist; - d_t_SFscaling d_SFscaling; - t_SFscaling SFscaling; - - // NN Kokkos::Views - d_t_bias bias; - d_t_weights weights; - t_bias h_bias; - t_weights h_weights; - int numLayers, numHiddenLayers, maxNeurons; - d_t_int numNeuronsPerLayer; - h_t_int h_numNeuronsPerLayer; - d_t_int AF; - h_t_int h_AF; - - // Per type Kokkos::Views - h_t_mass atomicEnergyOffset; - h_t_int h_numSFperElem; - d_t_int numSFperElem; - h_t_int h_numSFGperElem; - d_t_int numSFGperElem; - int maxSFperElem; - - using Mode::numElements; - using Mode::minNeighbors; - using Mode::minCutoffRadius; - using Mode::maxCutoffRadius; - using Mode::cutoffAlpha; - double meanEnergy; - using Mode::convEnergy; - using Mode::convLength; - ScalingType scalingType; - using Mode::settings; - using Mode::cutoffType; - std::vector elements; - std::vector elementStrings; }; -////////////////////////////////// -// Inlined function definitions // -////////////////////////////////// - -template -KOKKOS_INLINE_FUNCTION void -ModeCabana::compute_cutoff(CutoffFunction::CutoffType cutoffType, - double cutoffAlpha, double &fc, double &dfc, - double r, double rc, bool derivative) const -{ - double temp; - if (cutoffType == CutoffFunction::CT_TANHU) { - temp = tanh(1.0 - r / rc); - fc = temp * temp * temp; - if (derivative) - dfc = 3.0 * temp * temp * (temp * temp - 1.0) / rc; - } - - if (cutoffType == CutoffFunction::CT_COS) { - - double rci = rc * cutoffAlpha; - double iw = 1.0 / (rc - rci); - double PI = 4.0 * atan(1.0); - if (r < rci) { - fc = 1.0; - dfc = 0.0; - } else { - temp = cos(PI * (r - rci) * iw); - fc = 0.5 * (temp + 1.0); - if (derivative) - dfc = -0.5 * iw * PI * sqrt(1.0 - temp * temp); - } - } -} - -template -KOKKOS_INLINE_FUNCTION double -ModeCabana::scale(int attype, double value, int k, - d_t_SFscaling SFscaling_) const -{ - double scalingType = SFscaling_(attype, k, 7); - double scalingFactor = SFscaling_(attype, k, 6); - double Gmin = SFscaling_(attype, k, 0); - // double Gmax = SFscaling_(attype,k,1); - double Gmean = SFscaling_(attype, k, 2); - // double Gsigma = SFscaling_(attype,k,3); - double Smin = SFscaling_(attype, k, 4); - // double Smax = SFscaling_(attype,k,5); - - if (scalingType == 0.0) { - return value; - } else if (scalingType == 1.0) { - return Smin + scalingFactor * (value - Gmin); - } else if (scalingType == 2.0) { - return value - Gmean; - } else if (scalingType == 3.0) { - return Smin + scalingFactor * (value - Gmean); - } else if (scalingType == 4.0) { - return Smin + scalingFactor * (value - Gmean); - } else { - return 0.0; - } -} - } #include "ModeCabana_impl.h" diff --git a/src/libnnpif/CabanaMD/ModeCabana_impl.h b/src/libnnpif/CabanaMD/ModeCabana_impl.h index 068ab04ce..87b818af9 100644 --- a/src/libnnpif/CabanaMD/ModeCabana_impl.h +++ b/src/libnnpif/CabanaMD/ModeCabana_impl.h @@ -29,578 +29,6 @@ using namespace std; namespace nnp { -// TODO: call base, then copy to Views -template -void ModeCabana::setupElementMap() -{ - log << "\n"; - log << "*** SETUP: ELEMENT MAP ******************" - "**************************************\n"; - log << "\n"; - - elementStrings = split( reduce( settings["elements"] ) ); - - log << strpr( "Number of element strings found: %d\n", - elementStrings.size() ); - for ( size_t i = 0; i < elementStrings.size(); ++i ) - { - log << strpr( "Element %2zu: %2s\n", i, - elementStrings[i].c_str() ); - } - // resize to match number of element types - numElements = elementStrings.size(); - - log << "*****************************************" - "**************************************\n"; - - return; -} - -template -void ModeCabana::setupElements() -{ - log << "\n"; - log << "*** SETUP: ELEMENTS *********************" - "**************************************\n"; - log << "\n"; - - numElements = (size_t)atoi( settings["number_of_elements"].c_str() ); - atomicEnergyOffset = - h_t_mass( "Mode::atomicEnergyOffset", numElements ); - if ( numElements != elementStrings.size() ) - { - throw runtime_error( "ERROR: Inconsistent number of elements.\n" ); - } - log << strpr( "Number of elements is consistent: %zu\n", numElements ); - - for ( size_t i = 0; i < numElements; ++i ) - { - elements.push_back( ElementCabana( i ) ); - } - - if ( settings.keywordExists( "atom_energy" ) ) - { - Settings::KeyRange r = settings.getValues( "atom_energy" ); - for ( Settings::KeyMap::const_iterator it = r.first; - it != r.second; ++it ) - { - vector args = split( reduce( it->second.first ) ); - const char *estring = args.at( 0 ).c_str(); - for ( size_t i = 0; i < elementStrings.size(); ++i ) - { - if ( strcmp( elementStrings[i].c_str(), estring ) == 0 ) - atomicEnergyOffset( i ) = atof( args.at( 1 ).c_str() ); - } - } - } - - log << "Atomic energy offsets per element:\n"; - for ( size_t i = 0; i < elementStrings.size(); ++i ) - { - log << strpr( "Element %2zu: %16.8E\n", i, - atomicEnergyOffset( i ) ); - } - - log << "Energy offsets are automatically subtracted from reference " - "energies.\n"; - log << "*****************************************" - "**************************************\n"; - - return; -} - -template -void ModeCabana::setupSymmetryFunctions() -{ - maxSFperElem = 0; - h_numSFperElem = - h_t_int( "Mode::numSymmetryFunctionsPerElement", numElements ); - log << "\n"; - log << "*** SETUP: SYMMETRY FUNCTIONS ***********" - "**************************************\n"; - log << "\n"; - - // Only count SF per element; parse and add later - Settings::KeyRange r = settings.getValues( "symfunction_short" ); - for ( Settings::KeyMap::const_iterator it = r.first; it != r.second; - ++it ) - { - vector args = split( reduce( it->second.first ) ); - int type = 0; - const char *estring = args.at( 0 ).c_str(); - for ( size_t i = 0; i < elementStrings.size(); ++i ) - { - if ( strcmp( elementStrings[i].c_str(), estring ) == 0 ) - type = i; - } - h_numSFperElem( type )++; - - if ( h_numSFperElem( type ) > maxSFperElem ) - maxSFperElem = h_numSFperElem( type ); - } - Kokkos::deep_copy( h_numSFperElem, 0 ); - - // setup SF host views - // create device mirrors if needed below - SF = t_SF( "SymmetryFunctions", numElements, maxSFperElem ); - SFscaling = t_SFscaling( "SFscaling", numElements, maxSFperElem ); - // +1 to store size of memberlist - SFGmemberlist = t_SFGmemberlist( "SFGmemberlist", numElements, - maxSFperElem + 1, maxSFperElem + 1 ); - - r = settings.getValues( "symfunction_short" ); - for ( Settings::KeyMap::const_iterator it = r.first; it != r.second; - ++it ) - { - vector args = split( reduce( it->second.first ) ); - int type = 0; - const char *estring = args.at( 0 ).c_str(); - for ( size_t i = 0; i < elementStrings.size(); ++i ) - { - if ( strcmp( elementStrings[i].c_str(), estring ) == 0 ) - type = i; - } - elements.at( type ).addSymmetryFunction( it->second.first, - elementStrings, type, SF, - convLength, h_numSFperElem ); - } - - log << "Abbreviations:\n"; - log << "--------------\n"; - log << "ind .... Symmetry function index.\n"; - log << "ec ..... Central atom element.\n"; - log << "ty ..... Symmetry function type.\n"; - log << "e1 ..... Neighbor 1 element.\n"; - log << "e2 ..... Neighbor 2 element.\n"; - log << "eta .... Gaussian width eta.\n"; - log << "rs ..... Shift distance of Gaussian.\n"; - log << "la ..... Angle prefactor lambda.\n"; - log << "zeta ... Angle term exponent zeta.\n"; - log << "rc ..... Cutoff radius.\n"; - log << "ct ..... Cutoff type.\n"; - log << "ca ..... Cutoff alpha.\n"; - log << "ln ..... Line number in settings file.\n"; - log << "\n"; - maxCutoffRadius = 0.0; - - for ( vector::iterator it = elements.begin(); it != elements.end(); - ++it ) - { - int attype = it->getIndex(); - it->sortSymmetryFunctions( SF, h_numSFperElem, attype ); - maxCutoffRadius = - max( it->getMaxCutoffRadius( SF, attype, h_numSFperElem ), - maxCutoffRadius ); - it->setCutoffFunction( cutoffType, cutoffAlpha, SF, attype, - h_numSFperElem ); - log << strpr( - "Short range atomic symmetry functions element %2s :\n", - it->getSymbol().c_str() ); - log << "-----------------------------------------" - "--------------------------------------\n"; - log << " ind ec ty e1 e2 eta rs la " - "zeta rc ct ca ln\n"; - log << "-----------------------------------------" - "--------------------------------------\n"; - log << it->infoSymmetryFunctionParameters( SF, attype, h_numSFperElem ); - log << "-----------------------------------------" - "--------------------------------------\n"; - } - minNeighbors.resize( numElements, 0 ); - minCutoffRadius.resize( numElements, maxCutoffRadius ); - for ( size_t i = 0; i < numElements; ++i ) - { - int attype = elements.at( i ).getIndex(); - int nSF = h_numSFperElem( attype ); - minNeighbors.at( i ) = - elements.at( i ).getMinNeighbors( attype, SF, nSF ); - minCutoffRadius.at( i ) = - elements.at( i ).getMinCutoffRadius( SF, attype, h_numSFperElem ); - log << strpr( "Minimum cutoff radius for element %2s: %f\n", - elements.at( i ).getSymbol().c_str(), - minCutoffRadius.at( i ) / convLength ); - } - log << strpr( "Maximum cutoff radius (global) : %f\n", - maxCutoffRadius / convLength ); - - log << "*****************************************" - "**************************************\n"; - - numSFperElem = - Kokkos::create_mirror_view_and_copy( memory_space(), h_numSFperElem ); - - return; -} - -// TODO: call base, then copy to View -template -void ModeCabana::setupSymmetryFunctionScaling( string const &fileName ) -{ - log << "\n"; - log << "*** SETUP: SYMMETRY FUNCTION SCALING ****" - "**************************************\n"; - log << "\n"; - - log << "Equal scaling type for all symmetry functions:\n"; - if ( ( settings.keywordExists( "scale_symmetry_functions" ) ) && - ( !settings.keywordExists( "center_symmetry_functions" ) ) ) - { - scalingType = ST_SCALE; - log << strpr( "Scaling type::ST_SCALE (%d)\n", scalingType ); - log << "Gs = Smin + (Smax - Smin) * (G - Gmin) / (Gmax - Gmin)\n"; - } - else if ( ( !settings.keywordExists( "scale_symmetry_functions" ) ) && - ( settings.keywordExists( "center_symmetry_functions" ) ) ) - { - scalingType = ST_CENTER; - log << strpr( "Scaling type::ST_CENTER (%d)\n", scalingType ); - log << "Gs = G - Gmean\n"; - } - else if ( ( settings.keywordExists( "scale_symmetry_functions" ) ) && - ( settings.keywordExists( "center_symmetry_functions" ) ) ) - { - scalingType = ST_SCALECENTER; - log << strpr( "Scaling type::ST_SCALECENTER (%d)\n", scalingType ); - log << "Gs = Smin + (Smax - Smin) * (G - Gmean) / (Gmax - Gmin)\n"; - } - else if ( settings.keywordExists( "scale_symmetry_functions_sigma" ) ) - { - scalingType = ST_SCALESIGMA; - log << strpr( "Scaling type::ST_SCALESIGMA (%d)\n", scalingType ); - log << "Gs = Smin + (Smax - Smin) * (G - Gmean) / Gsigma\n"; - } - else - { - scalingType = ST_NONE; - log << strpr( "Scaling type::ST_NONE (%d)\n", scalingType ); - log << "Gs = G\n"; - log << "WARNING: No symmetry function scaling!\n"; - } - - double Smin = 0.0; - double Smax = 0.0; - if ( scalingType == ST_SCALE || scalingType == ST_SCALECENTER || - scalingType == ST_SCALESIGMA ) - { - if ( settings.keywordExists( "scale_min_short" ) ) - { - Smin = atof( settings["scale_min_short"].c_str() ); - } - else - { - log << "WARNING: Keyword \"scale_min_short\" not found.\n"; - log << " Default value for Smin = 0.0.\n"; - Smin = 0.0; - } - - if ( settings.keywordExists( "scale_max_short" ) ) - { - Smax = atof( settings["scale_max_short"].c_str() ); - } - else - { - log << "WARNING: Keyword \"scale_max_short\" not found.\n"; - log << " Default value for Smax = 1.0.\n"; - Smax = 1.0; - } - - log << strpr( "Smin = %f\n", Smin ); - log << strpr( "Smax = %f\n", Smax ); - } - - log << strpr( "Symmetry function scaling statistics from file: %s\n", - fileName.c_str() ); - log << "-----------------------------------------" - "--------------------------------------\n"; - ifstream file; - file.open( fileName.c_str() ); - if ( !file.is_open() ) - { - throw runtime_error( "ERROR: Could not open file: \"" + fileName + - "\".\n" ); - } - string line; - vector lines; - while ( getline( file, line ) ) - { - if ( line.at( 0 ) != '#' ) - lines.push_back( line ); - } - file.close(); - - log << "\n"; - log << "Abbreviations:\n"; - log << "--------------\n"; - log << "ind ..... Symmetry function index.\n"; - log << "min ..... Minimum symmetry function value.\n"; - log << "max ..... Maximum symmetry function value.\n"; - log << "mean .... Mean symmetry function value.\n"; - log << "sigma ... Standard deviation of symmetry function values.\n"; - log << "sf ...... Scaling factor for derivatives.\n"; - log << "Smin .... Desired minimum scaled symmetry function value.\n"; - 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(); - ++it ) - { - int attype = it->getIndex(); - it->setScaling( scalingType, lines, Smin, Smax, SF, SFscaling, attype, - h_numSFperElem ); - log << strpr( - "Scaling data for symmetry functions element %2s :\n", - it->getSymbol().c_str() ); - log << "-----------------------------------------" - "--------------------------------------\n"; - log << " ind min max mean sigma sf Smin " - "Smax t\n"; - log << "-----------------------------------------" - "--------------------------------------\n"; - log << it->infoSymmetryFunctionScaling( scalingType, SF, SFscaling, - attype, h_numSFperElem ); - log << "-----------------------------------------" - "--------------------------------------\n"; - lines.erase( lines.begin(), - lines.begin() + - it->numSymmetryFunctions( attype, h_numSFperElem ) ); - } - - log << "*****************************************" - "**************************************\n"; - - d_SF = Kokkos::create_mirror_view_and_copy( memory_space(), SF ); - d_SFscaling = - Kokkos::create_mirror_view_and_copy( memory_space(), SFscaling ); - d_SFGmemberlist = - Kokkos::create_mirror_view_and_copy( memory_space(), SFGmemberlist ); - - return; -} - -template -void ModeCabana::setupSymmetryFunctionGroups() -{ - log << "\n"; - log << "*** SETUP: SYMMETRY FUNCTION GROUPS *****" - "**************************************\n"; - log << "\n"; - - log << "Abbreviations:\n"; - log << "--------------\n"; - log << "ind .... Symmetry function group index.\n"; - log << "ec ..... Central atom element.\n"; - log << "ty ..... Symmetry function type.\n"; - log << "e1 ..... Neighbor 1 element.\n"; - log << "e2 ..... Neighbor 2 element.\n"; - log << "eta .... Gaussian width eta.\n"; - log << "rs ..... Shift distance of Gaussian.\n"; - log << "la ..... Angle prefactor lambda.\n"; - log << "zeta ... Angle term exponent zeta.\n"; - log << "rc ..... Cutoff radius.\n"; - log << "ct ..... Cutoff type.\n"; - log << "ca ..... Cutoff alpha.\n"; - log << "ln ..... Line number in settings file.\n"; - log << "mi ..... Member index.\n"; - log << "sfi .... Symmetry function index.\n"; - log << "e ...... Recalculate exponential term.\n"; - log << "\n"; - - h_numSFGperElem = - h_t_int( "Mode::numSymmetryFunctionGroupsPerElement", numElements ); - - for ( vector::iterator it = elements.begin(); it != elements.end(); - ++it ) - { - int attype = it->getIndex(); - it->setupSymmetryFunctionGroups( SF, SFGmemberlist, attype, - h_numSFperElem, h_numSFGperElem, - maxSFperElem ); - log << strpr( "Short range atomic symmetry function groups " - "element %2s :\n", - it->getSymbol().c_str() ); - log << "-----------------------------------------" - "--------------------------------------\n"; - log << " ind ec ty e1 e2 eta rs la " - "zeta rc ct ca ln mi sfi e\n"; - log << "-----------------------------------------" - "--------------------------------------\n"; - log << it->infoSymmetryFunctionGroups( SF, SFGmemberlist, attype, - h_numSFGperElem ); - log << "-----------------------------------------" - "--------------------------------------\n"; - } - - log << "*****************************************" - "**************************************\n"; - - numSFGperElem = - Kokkos::create_mirror_view_and_copy( memory_space(), h_numSFGperElem ); - - return; -} - -template -void ModeCabana::setupNeuralNetwork() -{ - log << "\n"; - log << "*** SETUP: NEURAL NETWORKS **************" - "**************************************\n"; - log << "\n"; - - numLayers = 2 + atoi( settings["global_hidden_layers_short"].c_str() ); - numHiddenLayers = numLayers - 2; - - h_numNeuronsPerLayer = h_t_int( "Mode::numNeuronsPerLayer", numLayers ); - h_AF = h_t_int( "Mode::ActivationFunctions", numLayers ); - - vector numNeuronsPerHiddenLayer = - split( reduce( settings["global_nodes_short"] ) ); - vector activationFunctions = - split( reduce( settings["global_activation_short"] ) ); - - for ( int i = 0; i < numLayers; i++ ) - { - if ( i == 0 ) - h_AF( i ) = 0; - else if ( i == numLayers - 1 ) - { - h_numNeuronsPerLayer( i ) = 1; - h_AF( i ) = 0; - } - else - { - h_numNeuronsPerLayer( i ) = - atoi( numNeuronsPerHiddenLayer.at( i - 1 ).c_str() ); - h_AF( i ) = 1; // TODO: hardcoded atoi(activationFunctions.at(i-1)); - } - } - - // TODO: add normalization of neurons - bool normalizeNeurons = settings.keywordExists( "normalize_nodes" ); - log << strpr( "Normalize neurons (all elements): %d\n", - (int)normalizeNeurons ); - log << "-----------------------------------------" - "--------------------------------------\n"; - - for ( vector::iterator it = elements.begin(); it != elements.end(); - ++it ) - { - int attype = it->getIndex(); - h_numNeuronsPerLayer( 0 ) = - it->numSymmetryFunctions( attype, h_numSFperElem ); - log << strpr( "Atomic short range NN for " - "element %2s :\n", - it->getSymbol().c_str() ); - - int numWeights = 0, numBiases = 0, numConnections = 0; - for ( int j = 1; j < numLayers; ++j ) - { - numWeights += - h_numNeuronsPerLayer( j - 1 ) * h_numNeuronsPerLayer( j ); - numBiases += h_numNeuronsPerLayer( j ); - } - numConnections = numWeights + numBiases; - log << strpr( "Number of weights : %6zu\n", numWeights ); - log << strpr( "Number of biases : %6zu\n", numBiases ); - log << strpr( "Number of connections: %6zu\n", numConnections ); - log << strpr( "Architecture " ); - for ( int j = 0; j < numLayers; ++j ) - log << strpr( " %4d", h_numNeuronsPerLayer( j ) ); - - log << "\n-----------------------------------------" - "--------------------------------------\n"; - } - - // initialize Views - maxNeurons = 0; - for ( int j = 0; j < numLayers; ++j ) - maxNeurons = max( maxNeurons, h_numNeuronsPerLayer( j ) ); - - h_bias = t_bias( "Mode::biases", numElements, numLayers, maxNeurons ); - h_weights = t_weights( "Mode::weights", numElements, numLayers, - maxNeurons, maxNeurons ); - - log << "*****************************************" - "**************************************\n"; - - return; -} - -template -void ModeCabana::setupNeuralNetworkWeights( string const &fileNameFormat ) -{ - log << "\n"; - log << "*** SETUP: NEURAL NETWORK WEIGHTS *******" - "**************************************\n"; - log << "\n"; - - log << strpr( "Weight file name format: %s\n", - fileNameFormat.c_str() ); - int count = 0; - int AN = 0; - for ( vector::iterator it = elements.begin(); it != elements.end(); - ++it ) - { - const char *estring = elementStrings[count].c_str(); - for ( size_t i = 0; i < knownElements.size(); ++i ) - { - if ( strcmp( knownElements[i].c_str(), estring ) == 0 ) - { - AN = i + 1; - break; - } - } - - string fileName = strpr( fileNameFormat.c_str(), AN ); - log << strpr( "Weight file for element %2s: %s\n", - elementStrings[count].c_str(), fileName.c_str() ); - ifstream file; - file.open( fileName.c_str() ); - if ( !file.is_open() ) - { - throw runtime_error( "ERROR: Could not open file: \"" + fileName + - "\".\n" ); - } - string line; - int attype = it->getIndex(); - int layer, start, end; - while ( getline( file, line ) ) - { - if ( line.at( 0 ) != '#' ) - { - vector splitLine = split( reduce( line ) ); - if ( strcmp( splitLine.at( 1 ).c_str(), "a" ) == 0 ) - { - layer = atoi( splitLine.at( 3 ).c_str() ); - start = atoi( splitLine.at( 4 ).c_str() ) - 1; - end = atoi( splitLine.at( 6 ).c_str() ) - 1; - h_weights( attype, layer, end, start ) = - atof( splitLine.at( 0 ).c_str() ); - } - else if ( strcmp( splitLine.at( 1 ).c_str(), "b" ) == 0 ) - { - layer = atoi( splitLine.at( 3 ).c_str() ) - 1; - start = atoi( splitLine.at( 4 ).c_str() ) - 1; - h_bias( attype, layer, start ) = - atof( splitLine.at( 0 ).c_str() ); - } - } - } - file.close(); - count += 1; - } - log << "*****************************************" - "**************************************\n"; - - bias = Kokkos::create_mirror_view_and_copy( memory_space(), h_bias ); - weights = Kokkos::create_mirror_view_and_copy( memory_space(), h_weights ); - AF = Kokkos::create_mirror_view_and_copy( memory_space(), h_AF ); - numNeuronsPerLayer = Kokkos::create_mirror_view_and_copy( - memory_space(), h_numNeuronsPerLayer ); - - return; -} - template template @@ -653,8 +81,8 @@ void ModeCabana::calculateSymmetryFunctionGroups( rij = sqrt( r2ij ); if ( e1 == nej && rij < rc ) { - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, false ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, false ); for ( size_t k = 0; k < size; ++k ) { memberindex = SFGmemberlist_( attype, groupIndex, k ); @@ -708,8 +136,8 @@ void ModeCabana::calculateSymmetryFunctionGroups( if ( ( e1 == nej || e2 == nej ) && rij < rc ) { // Calculate cutoff function and derivative. - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, false ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, false ); nek = type( k ); @@ -733,12 +161,12 @@ void ModeCabana::calculateSymmetryFunctionGroups( if ( r2jk < rc * rc ) { // Energy calculation. - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcik, pdfcik, rik, rc, false ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcik, pdfcik, rik, rc, false ); rjk = sqrt( r2jk ); - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcjk, pdfcjk, rjk, rc, false ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcjk, pdfcjk, rjk, rc, false ); double const rinvijik = 1.0 / rij / rik; double const costijk = @@ -821,7 +249,7 @@ void ModeCabana::calculateSymmetryFunctionGroups( pow( 2, ( 1 - SF_( attype, memberindex, 6 ) ) ); G( i, globalIndex ) = - scale( attype, raw_value, memberindex, SFscaling_ ); + this->scale( attype, raw_value, memberindex, SFscaling_ ); } } }; @@ -830,101 +258,6 @@ void ModeCabana::calculateSymmetryFunctionGroups( Kokkos::fence(); } -template -template -void ModeCabana::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 ); - auto dfdx = d_t_NN( "Mode::dfdx", N_local, numLayers, maxNeurons ); - auto inner = d_t_NN( "Mode::inner", N_local, numHiddenLayers, maxNeurons ); - auto outer = d_t_NN( "Mode::outer", N_local, numHiddenLayers, maxNeurons ); - - Kokkos::RangePolicy policy( 0, N_local ); - - // Create local copies for lambda - auto numSFperElem_ = numSFperElem; - auto numNeuronsPerLayer_ = numNeuronsPerLayer; - auto numLayers_ = numLayers; - auto numHiddenLayers_ = numHiddenLayers; - auto AF_ = AF; - auto weights_ = weights; - auto bias_ = bias; - - auto calc_nn_op = KOKKOS_LAMBDA( const int atomindex ) - { - int attype = type( atomindex ); - // set input layer of NN - int layer_0, layer_lminusone; - layer_0 = (int)numSFperElem_( attype ); - - for ( int k = 0; k < layer_0; ++k ) - NN( atomindex, 0, k ) = G( atomindex, k ); - // forward propagation - for ( int l = 1; l < numLayers_; l++ ) - { - if ( l == 1 ) - layer_lminusone = layer_0; - else - layer_lminusone = numNeuronsPerLayer_( l - 1 ); - double dtmp; - for ( int i = 0; i < numNeuronsPerLayer_( l ); i++ ) - { - dtmp = 0.0; - for ( int j = 0; j < layer_lminusone; j++ ) - dtmp += weights_( attype, l - 1, i, j ) * - NN( atomindex, l - 1, j ); - dtmp += bias_( attype, l - 1, i ); - if ( AF_( l ) == 0 ) - { - NN( atomindex, l, i ) = dtmp; - dfdx( atomindex, l, i ) = 1.0; - } - else if ( AF_( l ) == 1 ) - { - dtmp = tanh( dtmp ); - NN( atomindex, l, i ) = dtmp; - dfdx( atomindex, l, i ) = 1.0 - dtmp * dtmp; - } - } - } - - E( atomindex ) = NN( atomindex, numLayers_ - 1, 0 ); - - // derivative of network w.r.t NN inputs - for ( int k = 0; k < numNeuronsPerLayer_( 0 ); k++ ) - { - for ( int i = 0; i < numNeuronsPerLayer_( 1 ); i++ ) - inner( atomindex, 0, i ) = - weights_( attype, 0, i, k ) * dfdx( atomindex, 1, i ); - - for ( int l = 1; l < numHiddenLayers_ + 1; l++ ) - { - for ( int i2 = 0; i2 < numNeuronsPerLayer_( l + 1 ); i2++ ) - { - outer( atomindex, l - 1, i2 ) = 0.0; - - for ( int i1 = 0; i1 < numNeuronsPerLayer_( l ); i1++ ) - outer( atomindex, l - 1, i2 ) += - weights_( attype, l, i2, i1 ) * - inner( atomindex, l - 1, i1 ); - outer( atomindex, l - 1, i2 ) *= - dfdx( atomindex, l + 1, i2 ); - - if ( l < numHiddenLayers_ ) - inner( atomindex, l, i2 ) = - outer( atomindex, l - 1, i2 ); - } - } - dEdG( atomindex, k ) = outer( atomindex, numHiddenLayers_ - 1, 0 ); - } - }; - Kokkos::parallel_for( "Mode::calculateAtomicNeuralNetworks", policy, - calc_nn_op ); - Kokkos::fence(); -} - template template ::calculateForces( { // Energy calculation. // Calculate cutoff function and derivative. - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, true ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, true ); for ( size_t k = 0; k < size; ++k ) { globalIndex = SF_( @@ -1055,8 +388,8 @@ void ModeCabana::calculateForces( if ( ( e1 == nej || e2 == nej ) && rij < rc ) { // Calculate cutoff function and derivative. - compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, - rij, rc, true ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, pfcij, pdfcij, + rij, rc, true ); nek = type( k ); if ( ( e1 == nej && e2 == nek ) || @@ -1079,12 +412,12 @@ void ModeCabana::calculateForces( if ( r2jk < rc * rc ) { // Energy calculation. - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcik, pdfcik, rik, rc, true ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcik, pdfcik, rik, rc, true ); rjk = sqrt( r2jk ); - compute_cutoff( cutoffType_, cutoffAlpha_, - pfcjk, pdfcjk, rjk, rc, true ); + this->compute_cutoff( cutoffType_, cutoffAlpha_, + pfcjk, pdfcjk, rjk, rc, true ); double const rinvijik = 1.0 / rij / rik; double const costijk = diff --git a/src/libnnpif/CabanaMD/ElementCabana.h b/src/libnnpif/Kokkos/ElementKokkos.h similarity index 96% rename from src/libnnpif/CabanaMD/ElementCabana.h rename to src/libnnpif/Kokkos/ElementKokkos.h index 69328dd8d..09180dece 100644 --- a/src/libnnpif/CabanaMD/ElementCabana.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/CabanaMD/ElementCabana_impl.h b/src/libnnpif/Kokkos/ElementKokkos_impl.h similarity index 94% rename from src/libnnpif/CabanaMD/ElementCabana_impl.h rename to src/libnnpif/Kokkos/ElementKokkos_impl.h index 3feae0507..f929e9417 100644 --- a/src/libnnpif/CabanaMD/ElementCabana_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,18 +31,17 @@ using namespace std; namespace nnp { -ElementCabana::ElementCabana( size_t const _index ) +inline ElementKokkos::ElementKokkos( size_t const _index ) : Element() { index = _index; atomicEnergyOffset = 0.0; - neuralNetwork = NULL; } -ElementCabana::~ElementCabana() {} +inline 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 +168,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 +204,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 +267,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 +292,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 +307,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 +389,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 +411,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 +424,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 +444,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 +464,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 +476,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 +489,9 @@ double ElementCabana::getMaxCutoffRadius( t_SF SF, int attype, // TODO:add functionality /* -void ElementCabana::updateSymmetryFunctionStatistics +void ElementKokkos::updateSymmetryFunctionStatistics */ } + +#endif diff --git a/src/libnnpif/Kokkos/ModeKokkos.h b/src/libnnpif/Kokkos/ModeKokkos.h new file mode 100644 index 000000000..0708b084c --- /dev/null +++ b/src/libnnpif/Kokkos/ModeKokkos.h @@ -0,0 +1,340 @@ +// 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_KOKKOS_H +#define MODE_KOKKOS_H + +#include "ElementKokkos.h" +#include "typesKokkos.h" + +#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 Kokkos main NNP class. + * + * The main n2p2 functions for computing energies and forces are replaced + * to use the Kokkos library. Most setup functions are + * overridden; some are replaced as needed to work within device kernels. + */ +template +class ModeKokkos : public Mode +{ + + public: + using Mode::Mode; + + // Kokkos settings + using device_type = t_device; + using exe_space = typename device_type::execution_space; + using memory_space = typename device_type::memory_space; + typedef typename exe_space::array_layout layout; + using host_space = Kokkos::HostSpace; + + // Per-type Kokkos::Views + using d_t_mass = Kokkos::View; + using h_t_mass = Kokkos::View; + using d_t_int = Kokkos::View; + using h_t_int = Kokkos::View; + + // SymmetryFunctionTypes Kokkos::Views + using d_t_SF = Kokkos::View; + using t_SF = Kokkos::View; + using d_t_SFscaling = Kokkos::View; + using t_SFscaling = Kokkos::View; + using d_t_SFGmemberlist = Kokkos::View; + using t_SFGmemberlist = Kokkos::View; + + // NN Kokkos::Views + using d_t_bias = Kokkos::View; + using t_bias = Kokkos::View; + using d_t_weights = Kokkos::View; + using t_weights = Kokkos::View; + using d_t_NN = Kokkos::View; + + /** Set up the element map. + * + * Uses keyword `elements`. This function should follow immediately after + * settings are loaded via loadSettingsFile(). + */ + void setupElementMap() override; + /** Set up all Element instances. + * + * Uses keywords `number_of_elements` and `atom_energy`. This function + * should follow immediately after setupElementMap(). + */ + void setupElements() override; + /** Set up all symmetry functions. + * + * Uses keyword `symfunction_short`. Reads all symmetry functions from + * settings and automatically assigns them to the correct element. + */ + void setupSymmetryFunctions() override; + /** Set up symmetry function scaling from file. + * + * @param[in] fileName Scaling file name. + * + * Uses keywords `scale_symmetry_functions`, `center_symmetry_functions`, + * `scale_symmetry_functions_sigma`, `scale_min_short` and + * `scale_max_short`. Reads in scaling information and sets correct scaling + * behavior for all symmetry functions. Call after + * setupSymmetryFunctions(). + */ + void setupSymmetryFunctionScaling( + std::string const &fileName = "scaling.data") override; + /** Set up symmetry function groups. + * + * Does not use any keywords. Call after setupSymmetryFunctions() and + * ensure that correct scaling behavior has already been set. + */ + void setupSymmetryFunctionGroups() override; + /** Set up neural networks for all elements. + * + * Uses keywords `global_hidden_layers_short`, `global_nodes_short`, + * `global_activation_short`, `normalize_nodes`. Call after + * setupSymmetryFunctions(), only then the number of input layer neurons is + * known. + */ + void setupNeuralNetwork() override; + /** Set up neural network weights from files. + * + * @param[in] fileNameFormat Format for weights file name. The string must + * contain one placeholder for the atomic number. + * + * Does not use any keywords. The weight files should contain one weight + * per line, see NeuralNetwork::setConnections() for the correct order. + */ + void setupNeuralNetworkWeights( + std::string const &fileNameFormatShort = "weights.%03zu.data", + std::string const &fileNameFormatCharge = "weightse.%03zu.data" ) override; + + /* Compute cutoff for a single atom pair. + * + * @param[in] cutoffType Cutoff type. + * @param[in] fc Cutoff function value. + * @param[in] dfc Derivative of cutoff function. + * @param[in] r Atom separation distance. + * @param[in] rc Cutoff radius. + * @param[in] derivative If the cutoff derivative should be computed. + * + * Callable within device kernel. + */ + KOKKOS_INLINE_FUNCTION + void compute_cutoff(CutoffFunction::CutoffType cutoffType, double cutoffAlpha, + double &fc, double &dfc, double r, double rc, + bool derivative) const; + /* + * Scale one symmetry function. + * + * @param[in] attype Atom type. + * @param[in] value Scaled symmetry function value. + * @param[in] k Symmetry function index. + * @param[in] SFscaling Kokkos host View of symmetry function scaling. + * + * Callable within device kernel. + */ + KOKKOS_INLINE_FUNCTION + double scale(int attype, double value, int k, d_t_SFscaling SFscaling) const; + + /** Calculate forces for all atoms in given structure. + * + * @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. + * + * 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 ); + + /** Calculate atomic neural networks for all atoms in given structure. + * + * @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 Kokkos View of energy per atom. + * @param[in] N_local Number of atoms. + * + * The atomic energy is stored in E. + */ + template + void calculateAtomicNeuralNetworks(t_slice_type type, t_slice_G G, + t_slice_dEdG dEdG, t_slice_E E, + int N_local); + + /** Calculate all symmetry function groups for all atoms in given + * structure. + * + * @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. + * + * 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); + + using Mode::log; + + /// list of element symbols in order of periodic table + // (duplicated from ElementMap) + std::vector knownElements = { + "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", + "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", + "Mn", "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", + "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", + "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", + "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", + "Ta", "W", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", + "At", "Rn", "Fr", "Ra", "Ac", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", + "Bk", "Cf", "Es", "Fm", "Md", "No"}; + + // Symmetry function Kokkos::Views + d_t_SF d_SF; + t_SF SF; + d_t_SFGmemberlist d_SFGmemberlist; + t_SFGmemberlist SFGmemberlist; + d_t_SFscaling d_SFscaling; + t_SFscaling SFscaling; + + // NN Kokkos::Views + d_t_bias bias; + d_t_weights weights; + t_bias h_bias; + t_weights h_weights; + int numLayers, numHiddenLayers, maxNeurons; + d_t_int numNeuronsPerLayer; + h_t_int h_numNeuronsPerLayer; + d_t_int AF; + h_t_int h_AF; + + // Per type Kokkos::Views + h_t_mass atomicEnergyOffset; + h_t_int h_numSFperElem; + d_t_int numSFperElem; + h_t_int h_numSFGperElem; + d_t_int numSFGperElem; + int maxSFperElem; + + using Mode::numElements; + using Mode::minNeighbors; + using Mode::minCutoffRadius; + using Mode::maxCutoffRadius; + using Mode::cutoffAlpha; + double meanEnergy; + using Mode::convEnergy; + using Mode::convLength; + ScalingType scalingType; + using Mode::settings; + using Mode::cutoffType; + std::vector elements; + std::vector elementStrings; +}; + +////////////////////////////////// +// Inlined function definitions // +////////////////////////////////// + +template +KOKKOS_INLINE_FUNCTION void +ModeKokkos::compute_cutoff(CutoffFunction::CutoffType cutoffType, + double cutoffAlpha, double &fc, double &dfc, + double r, double rc, bool derivative) const +{ + double temp; + if (cutoffType == CutoffFunction::CT_TANHU) { + temp = tanh(1.0 - r / rc); + fc = temp * temp * temp; + if (derivative) + dfc = 3.0 * temp * temp * (temp * temp - 1.0) / rc; + } + + if (cutoffType == CutoffFunction::CT_COS) { + + double rci = rc * cutoffAlpha; + double iw = 1.0 / (rc - rci); + double PI = 4.0 * atan(1.0); + if (r < rci) { + fc = 1.0; + dfc = 0.0; + } else { + temp = cos(PI * (r - rci) * iw); + fc = 0.5 * (temp + 1.0); + if (derivative) + dfc = -0.5 * iw * PI * sqrt(1.0 - temp * temp); + } + } +} + +template +KOKKOS_INLINE_FUNCTION double +ModeKokkos::scale(int attype, double value, int k, + d_t_SFscaling SFscaling_) const +{ + double scalingType = SFscaling_(attype, k, 7); + double scalingFactor = SFscaling_(attype, k, 6); + double Gmin = SFscaling_(attype, k, 0); + // double Gmax = SFscaling_(attype,k,1); + double Gmean = SFscaling_(attype, k, 2); + // double Gsigma = SFscaling_(attype,k,3); + double Smin = SFscaling_(attype, k, 4); + // double Smax = SFscaling_(attype,k,5); + + if (scalingType == 0.0) { + return value; + } else if (scalingType == 1.0) { + return Smin + scalingFactor * (value - Gmin); + } else if (scalingType == 2.0) { + return value - Gmean; + } else if (scalingType == 3.0) { + return Smin + scalingFactor * (value - Gmean); + } else if (scalingType == 4.0) { + return Smin + scalingFactor * (value - Gmean); + } else { + return 0.0; + } +} + +} + +#include "ModeKokkos_impl.h" + +#endif diff --git a/src/libnnpif/Kokkos/ModeKokkos_impl.h b/src/libnnpif/Kokkos/ModeKokkos_impl.h new file mode 100644 index 000000000..fd231d23c --- /dev/null +++ b/src/libnnpif/Kokkos/ModeKokkos_impl.h @@ -0,0 +1,1257 @@ +// 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_KOKKOS_IMPL_H +#define MODE_KOKKOS_IMPL_H + +#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 +{ + +// TODO: call base, then copy to Views +template +void ModeKokkos::setupElementMap() +{ + log << "\n"; + log << "*** SETUP: ELEMENT MAP ******************" + "**************************************\n"; + log << "\n"; + + elementStrings = split( reduce( settings["elements"] ) ); + + log << strpr( "Number of element strings found: %d\n", + elementStrings.size() ); + for ( size_t i = 0; i < elementStrings.size(); ++i ) + { + log << strpr( "Element %2zu: %2s\n", i, + elementStrings[i].c_str() ); + } + // resize to match number of element types + numElements = elementStrings.size(); + + log << "*****************************************" + "**************************************\n"; + + return; +} + +template +void ModeKokkos::setupElements() +{ + log << "\n"; + log << "*** SETUP: ELEMENTS *********************" + "**************************************\n"; + log << "\n"; + + numElements = (size_t)atoi( settings["number_of_elements"].c_str() ); + atomicEnergyOffset = + h_t_mass( "Mode::atomicEnergyOffset", numElements ); + if ( numElements != elementStrings.size() ) + { + throw runtime_error( "ERROR: Inconsistent number of elements.\n" ); + } + log << strpr( "Number of elements is consistent: %zu\n", numElements ); + + for ( size_t i = 0; i < numElements; ++i ) + { + elements.push_back( ElementKokkos( i ) ); + } + + if ( settings.keywordExists( "atom_energy" ) ) + { + Settings::KeyRange r = settings.getValues( "atom_energy" ); + for ( Settings::KeyMap::const_iterator it = r.first; + it != r.second; ++it ) + { + vector args = split( reduce( it->second.first ) ); + const char *estring = args.at( 0 ).c_str(); + for ( size_t i = 0; i < elementStrings.size(); ++i ) + { + if ( strcmp( elementStrings[i].c_str(), estring ) == 0 ) + atomicEnergyOffset( i ) = atof( args.at( 1 ).c_str() ); + } + } + } + + log << "Atomic energy offsets per element:\n"; + for ( size_t i = 0; i < elementStrings.size(); ++i ) + { + log << strpr( "Element %2zu: %16.8E\n", i, + atomicEnergyOffset( i ) ); + } + + log << "Energy offsets are automatically subtracted from reference " + "energies.\n"; + log << "*****************************************" + "**************************************\n"; + + return; +} + +template +void ModeKokkos::setupSymmetryFunctions() +{ + maxSFperElem = 0; + h_numSFperElem = + h_t_int( "Mode::numSymmetryFunctionsPerElement", numElements ); + log << "\n"; + log << "*** SETUP: SYMMETRY FUNCTIONS ***********" + "**************************************\n"; + log << "\n"; + + // Only count SF per element; parse and add later + Settings::KeyRange r = settings.getValues( "symfunction_short" ); + for ( Settings::KeyMap::const_iterator it = r.first; it != r.second; + ++it ) + { + vector args = split( reduce( it->second.first ) ); + int type = 0; + const char *estring = args.at( 0 ).c_str(); + for ( size_t i = 0; i < elementStrings.size(); ++i ) + { + if ( strcmp( elementStrings[i].c_str(), estring ) == 0 ) + type = i; + } + h_numSFperElem( type )++; + + if ( h_numSFperElem( type ) > maxSFperElem ) + maxSFperElem = h_numSFperElem( type ); + } + Kokkos::deep_copy( h_numSFperElem, 0 ); + + // setup SF host views + // create device mirrors if needed below + SF = t_SF( "SymmetryFunctions", numElements, maxSFperElem ); + SFscaling = t_SFscaling( "SFscaling", numElements, maxSFperElem ); + // +1 to store size of memberlist + SFGmemberlist = t_SFGmemberlist( "SFGmemberlist", numElements, + maxSFperElem + 1, maxSFperElem + 1 ); + + r = settings.getValues( "symfunction_short" ); + for ( Settings::KeyMap::const_iterator it = r.first; it != r.second; + ++it ) + { + vector args = split( reduce( it->second.first ) ); + int type = 0; + const char *estring = args.at( 0 ).c_str(); + for ( size_t i = 0; i < elementStrings.size(); ++i ) + { + if ( strcmp( elementStrings[i].c_str(), estring ) == 0 ) + type = i; + } + elements.at( type ).addSymmetryFunction( it->second.first, + elementStrings, type, SF, + convLength, h_numSFperElem ); + } + + log << "Abbreviations:\n"; + log << "--------------\n"; + log << "ind .... Symmetry function index.\n"; + log << "ec ..... Central atom element.\n"; + log << "ty ..... Symmetry function type.\n"; + log << "e1 ..... Neighbor 1 element.\n"; + log << "e2 ..... Neighbor 2 element.\n"; + log << "eta .... Gaussian width eta.\n"; + log << "rs ..... Shift distance of Gaussian.\n"; + log << "la ..... Angle prefactor lambda.\n"; + log << "zeta ... Angle term exponent zeta.\n"; + log << "rc ..... Cutoff radius.\n"; + log << "ct ..... Cutoff type.\n"; + log << "ca ..... Cutoff alpha.\n"; + log << "ln ..... Line number in settings file.\n"; + log << "\n"; + maxCutoffRadius = 0.0; + + for ( vector::iterator it = elements.begin(); it != elements.end(); + ++it ) + { + int attype = it->getIndex(); + it->sortSymmetryFunctions( SF, h_numSFperElem, attype ); + maxCutoffRadius = + max( it->getMaxCutoffRadius( SF, attype, h_numSFperElem ), + maxCutoffRadius ); + it->setCutoffFunction( cutoffType, cutoffAlpha, SF, attype, + h_numSFperElem ); + log << strpr( + "Short range atomic symmetry functions element %2s :\n", + it->getSymbol().c_str() ); + log << "-----------------------------------------" + "--------------------------------------\n"; + log << " ind ec ty e1 e2 eta rs la " + "zeta rc ct ca ln\n"; + log << "-----------------------------------------" + "--------------------------------------\n"; + log << it->infoSymmetryFunctionParameters( SF, attype, h_numSFperElem ); + log << "-----------------------------------------" + "--------------------------------------\n"; + } + minNeighbors.resize( numElements, 0 ); + minCutoffRadius.resize( numElements, maxCutoffRadius ); + for ( size_t i = 0; i < numElements; ++i ) + { + int attype = elements.at( i ).getIndex(); + int nSF = h_numSFperElem( attype ); + minNeighbors.at( i ) = + elements.at( i ).getMinNeighbors( attype, SF, nSF ); + minCutoffRadius.at( i ) = + elements.at( i ).getMinCutoffRadius( SF, attype, h_numSFperElem ); + log << strpr( "Minimum cutoff radius for element %2s: %f\n", + elements.at( i ).getSymbol().c_str(), + minCutoffRadius.at( i ) / convLength ); + } + log << strpr( "Maximum cutoff radius (global) : %f\n", + maxCutoffRadius / convLength ); + + log << "*****************************************" + "**************************************\n"; + + numSFperElem = + Kokkos::create_mirror_view_and_copy( memory_space(), h_numSFperElem ); + + return; +} + +// TODO: call base, then copy to View +template +void ModeKokkos::setupSymmetryFunctionScaling( string const &fileName ) +{ + log << "\n"; + log << "*** SETUP: SYMMETRY FUNCTION SCALING ****" + "**************************************\n"; + log << "\n"; + + log << "Equal scaling type for all symmetry functions:\n"; + if ( ( settings.keywordExists( "scale_symmetry_functions" ) ) && + ( !settings.keywordExists( "center_symmetry_functions" ) ) ) + { + scalingType = ST_SCALE; + log << strpr( "Scaling type::ST_SCALE (%d)\n", scalingType ); + log << "Gs = Smin + (Smax - Smin) * (G - Gmin) / (Gmax - Gmin)\n"; + } + else if ( ( !settings.keywordExists( "scale_symmetry_functions" ) ) && + ( settings.keywordExists( "center_symmetry_functions" ) ) ) + { + scalingType = ST_CENTER; + log << strpr( "Scaling type::ST_CENTER (%d)\n", scalingType ); + log << "Gs = G - Gmean\n"; + } + else if ( ( settings.keywordExists( "scale_symmetry_functions" ) ) && + ( settings.keywordExists( "center_symmetry_functions" ) ) ) + { + scalingType = ST_SCALECENTER; + log << strpr( "Scaling type::ST_SCALECENTER (%d)\n", scalingType ); + log << "Gs = Smin + (Smax - Smin) * (G - Gmean) / (Gmax - Gmin)\n"; + } + else if ( settings.keywordExists( "scale_symmetry_functions_sigma" ) ) + { + scalingType = ST_SCALESIGMA; + log << strpr( "Scaling type::ST_SCALESIGMA (%d)\n", scalingType ); + log << "Gs = Smin + (Smax - Smin) * (G - Gmean) / Gsigma\n"; + } + else + { + scalingType = ST_NONE; + log << strpr( "Scaling type::ST_NONE (%d)\n", scalingType ); + log << "Gs = G\n"; + log << "WARNING: No symmetry function scaling!\n"; + } + + double Smin = 0.0; + double Smax = 0.0; + if ( scalingType == ST_SCALE || scalingType == ST_SCALECENTER || + scalingType == ST_SCALESIGMA ) + { + if ( settings.keywordExists( "scale_min_short" ) ) + { + Smin = atof( settings["scale_min_short"].c_str() ); + } + else + { + log << "WARNING: Keyword \"scale_min_short\" not found.\n"; + log << " Default value for Smin = 0.0.\n"; + Smin = 0.0; + } + + if ( settings.keywordExists( "scale_max_short" ) ) + { + Smax = atof( settings["scale_max_short"].c_str() ); + } + else + { + log << "WARNING: Keyword \"scale_max_short\" not found.\n"; + log << " Default value for Smax = 1.0.\n"; + Smax = 1.0; + } + + log << strpr( "Smin = %f\n", Smin ); + log << strpr( "Smax = %f\n", Smax ); + } + + log << strpr( "Symmetry function scaling statistics from file: %s\n", + fileName.c_str() ); + log << "-----------------------------------------" + "--------------------------------------\n"; + ifstream file; + file.open( fileName.c_str() ); + if ( !file.is_open() ) + { + throw runtime_error( "ERROR: Could not open file: \"" + fileName + + "\".\n" ); + } + string line; + vector lines; + while ( getline( file, line ) ) + { + if ( line.at( 0 ) != '#' ) + lines.push_back( line ); + } + file.close(); + + log << "\n"; + log << "Abbreviations:\n"; + log << "--------------\n"; + log << "ind ..... Symmetry function index.\n"; + log << "min ..... Minimum symmetry function value.\n"; + log << "max ..... Maximum symmetry function value.\n"; + log << "mean .... Mean symmetry function value.\n"; + log << "sigma ... Standard deviation of symmetry function values.\n"; + log << "sf ...... Scaling factor for derivatives.\n"; + log << "Smin .... Desired minimum scaled symmetry function value.\n"; + 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(); + ++it ) + { + int attype = it->getIndex(); + it->setScaling( scalingType, lines, Smin, Smax, SF, SFscaling, attype, + h_numSFperElem ); + log << strpr( + "Scaling data for symmetry functions element %2s :\n", + it->getSymbol().c_str() ); + log << "-----------------------------------------" + "--------------------------------------\n"; + log << " ind min max mean sigma sf Smin " + "Smax t\n"; + log << "-----------------------------------------" + "--------------------------------------\n"; + log << it->infoSymmetryFunctionScaling( scalingType, SF, SFscaling, + attype, h_numSFperElem ); + log << "-----------------------------------------" + "--------------------------------------\n"; + lines.erase( lines.begin(), + lines.begin() + + it->numSymmetryFunctions( attype, h_numSFperElem ) ); + } + + log << "*****************************************" + "**************************************\n"; + + d_SF = Kokkos::create_mirror_view_and_copy( memory_space(), SF ); + d_SFscaling = + Kokkos::create_mirror_view_and_copy( memory_space(), SFscaling ); + d_SFGmemberlist = + Kokkos::create_mirror_view_and_copy( memory_space(), SFGmemberlist ); + + return; +} + +template +void ModeKokkos::setupSymmetryFunctionGroups() +{ + log << "\n"; + log << "*** SETUP: SYMMETRY FUNCTION GROUPS *****" + "**************************************\n"; + log << "\n"; + + log << "Abbreviations:\n"; + log << "--------------\n"; + log << "ind .... Symmetry function group index.\n"; + log << "ec ..... Central atom element.\n"; + log << "ty ..... Symmetry function type.\n"; + log << "e1 ..... Neighbor 1 element.\n"; + log << "e2 ..... Neighbor 2 element.\n"; + log << "eta .... Gaussian width eta.\n"; + log << "rs ..... Shift distance of Gaussian.\n"; + log << "la ..... Angle prefactor lambda.\n"; + log << "zeta ... Angle term exponent zeta.\n"; + log << "rc ..... Cutoff radius.\n"; + log << "ct ..... Cutoff type.\n"; + log << "ca ..... Cutoff alpha.\n"; + log << "ln ..... Line number in settings file.\n"; + log << "mi ..... Member index.\n"; + log << "sfi .... Symmetry function index.\n"; + log << "e ...... Recalculate exponential term.\n"; + log << "\n"; + + h_numSFGperElem = + h_t_int( "Mode::numSymmetryFunctionGroupsPerElement", numElements ); + + for ( vector::iterator it = elements.begin(); it != elements.end(); + ++it ) + { + int attype = it->getIndex(); + it->setupSymmetryFunctionGroups( SF, SFGmemberlist, attype, + h_numSFperElem, h_numSFGperElem, + maxSFperElem ); + log << strpr( "Short range atomic symmetry function groups " + "element %2s :\n", + it->getSymbol().c_str() ); + log << "-----------------------------------------" + "--------------------------------------\n"; + log << " ind ec ty e1 e2 eta rs la " + "zeta rc ct ca ln mi sfi e\n"; + log << "-----------------------------------------" + "--------------------------------------\n"; + log << it->infoSymmetryFunctionGroups( SF, SFGmemberlist, attype, + h_numSFGperElem ); + log << "-----------------------------------------" + "--------------------------------------\n"; + } + + log << "*****************************************" + "**************************************\n"; + + numSFGperElem = + Kokkos::create_mirror_view_and_copy( memory_space(), h_numSFGperElem ); + + return; +} + +template +void ModeKokkos::setupNeuralNetwork() +{ + log << "\n"; + log << "*** SETUP: NEURAL NETWORKS **************" + "**************************************\n"; + log << "\n"; + + numLayers = 2 + atoi( settings["global_hidden_layers_short"].c_str() ); + numHiddenLayers = numLayers - 2; + + h_numNeuronsPerLayer = h_t_int( "Mode::numNeuronsPerLayer", numLayers ); + h_AF = h_t_int( "Mode::ActivationFunctions", numLayers ); + + vector numNeuronsPerHiddenLayer = + split( reduce( settings["global_nodes_short"] ) ); + vector activationFunctions = + split( reduce( settings["global_activation_short"] ) ); + + for ( int i = 0; i < numLayers; i++ ) + { + if ( i == 0 ) + h_AF( i ) = 0; + else if ( i == numLayers - 1 ) + { + h_numNeuronsPerLayer( i ) = 1; + h_AF( i ) = 0; + } + else + { + h_numNeuronsPerLayer( i ) = + atoi( numNeuronsPerHiddenLayer.at( i - 1 ).c_str() ); + h_AF( i ) = 1; // TODO: hardcoded atoi(activationFunctions.at(i-1)); + } + } + + // TODO: add normalization of neurons + bool normalizeNeurons = settings.keywordExists( "normalize_nodes" ); + log << strpr( "Normalize neurons (all elements): %d\n", + (int)normalizeNeurons ); + log << "-----------------------------------------" + "--------------------------------------\n"; + + for ( vector::iterator it = elements.begin(); it != elements.end(); + ++it ) + { + int attype = it->getIndex(); + h_numNeuronsPerLayer( 0 ) = + it->numSymmetryFunctions( attype, h_numSFperElem ); + log << strpr( "Atomic short range NN for " + "element %2s :\n", + it->getSymbol().c_str() ); + + int numWeights = 0, numBiases = 0, numConnections = 0; + for ( int j = 1; j < numLayers; ++j ) + { + numWeights += + h_numNeuronsPerLayer( j - 1 ) * h_numNeuronsPerLayer( j ); + numBiases += h_numNeuronsPerLayer( j ); + } + numConnections = numWeights + numBiases; + log << strpr( "Number of weights : %6zu\n", numWeights ); + log << strpr( "Number of biases : %6zu\n", numBiases ); + log << strpr( "Number of connections: %6zu\n", numConnections ); + log << strpr( "Architecture " ); + for ( int j = 0; j < numLayers; ++j ) + log << strpr( " %4d", h_numNeuronsPerLayer( j ) ); + + log << "\n-----------------------------------------" + "--------------------------------------\n"; + } + + // initialize Views + maxNeurons = 0; + for ( int j = 0; j < numLayers; ++j ) + maxNeurons = max( maxNeurons, h_numNeuronsPerLayer( j ) ); + + h_bias = t_bias( "Mode::biases", numElements, numLayers, maxNeurons ); + h_weights = t_weights( "Mode::weights", numElements, numLayers, + maxNeurons, maxNeurons ); + + log << "*****************************************" + "**************************************\n"; + + return; +} + +template +void ModeKokkos::setupNeuralNetworkWeights( string const &fileNameFormat, + string const &) +{ + log << "\n"; + log << "*** SETUP: NEURAL NETWORK WEIGHTS *******" + "**************************************\n"; + log << "\n"; + + log << strpr( "Weight file name format: %s\n", + fileNameFormat.c_str() ); + int count = 0; + int AN = 0; + for ( vector::iterator it = elements.begin(); it != elements.end(); + ++it ) + { + const char *estring = elementStrings[count].c_str(); + for ( size_t i = 0; i < knownElements.size(); ++i ) + { + if ( strcmp( knownElements[i].c_str(), estring ) == 0 ) + { + AN = i + 1; + break; + } + } + + string fileName = strpr( fileNameFormat.c_str(), AN ); + log << strpr( "Weight file for element %2s: %s\n", + elementStrings[count].c_str(), fileName.c_str() ); + ifstream file; + file.open( fileName.c_str() ); + if ( !file.is_open() ) + { + throw runtime_error( "ERROR: Could not open file: \"" + fileName + + "\".\n" ); + } + string line; + int attype = it->getIndex(); + int layer, start, end; + while ( getline( file, line ) ) + { + if ( line.at( 0 ) != '#' ) + { + vector splitLine = split( reduce( line ) ); + if ( strcmp( splitLine.at( 1 ).c_str(), "a" ) == 0 ) + { + layer = atoi( splitLine.at( 3 ).c_str() ); + start = atoi( splitLine.at( 4 ).c_str() ) - 1; + end = atoi( splitLine.at( 6 ).c_str() ) - 1; + h_weights( attype, layer, end, start ) = + atof( splitLine.at( 0 ).c_str() ); + } + else if ( strcmp( splitLine.at( 1 ).c_str(), "b" ) == 0 ) + { + layer = atoi( splitLine.at( 3 ).c_str() ) - 1; + start = atoi( splitLine.at( 4 ).c_str() ) - 1; + h_bias( attype, layer, start ) = + atof( splitLine.at( 0 ).c_str() ); + } + } + } + file.close(); + count += 1; + } + log << "*****************************************" + "**************************************\n"; + + bias = Kokkos::create_mirror_view_and_copy( memory_space(), h_bias ); + weights = Kokkos::create_mirror_view_and_copy( memory_space(), h_weights ); + AF = Kokkos::create_mirror_view_and_copy( memory_space(), h_AF ); + numNeuronsPerLayer = Kokkos::create_mirror_view_and_copy( + memory_space(), h_numNeuronsPerLayer ); + + return; +} + +template +template +void ModeKokkos::calculateSymmetryFunctionGroups( + t_x x, t_type type, t_G G, t_neigh_list neigh_list, + int N_local ) +{ + Kokkos::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 ) + { + 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; + + 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 ) + { + 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_ ); + + for(int jj = 0; jj < num_neighs_i; jj++) + { + 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 ) + { + 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( "Mode::calculateRadialSymmetryFunctionGroups", + policy, calc_radial_symm_op ); + Kokkos::fence(); + + auto calc_angular_symm_op = KOKKOS_LAMBDA( const int i ) + { + 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; + + 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 ) + { + 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_ ); + + for(int jj = 0; jj < num_neighs_i; jj++) + { + 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 ) + { + // 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++) + { + T_INT k = neighs_i(kk); + 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. + 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( "Mode::calculateAngularSymmetryFunctionGroups", + policy, calc_angular_symm_op ); + 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 ) = + scale( attype, raw_value, memberindex, SFscaling_ ); + } + } + }; + Kokkos::parallel_for( "Mode::scaleSymmetryFunctionGroups", policy, + scale_symm_op ); + Kokkos::fence(); +} + +template +template +void ModeKokkos::calculateAtomicNeuralNetworks( + 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 ); + auto inner = d_t_NN( "Mode::inner", N_local, numHiddenLayers, maxNeurons ); + auto outer = d_t_NN( "Mode::outer", N_local, numHiddenLayers, maxNeurons ); + + Kokkos::RangePolicy policy( 0, N_local ); + + // Create local copies for lambda + auto numSFperElem_ = numSFperElem; + auto numNeuronsPerLayer_ = numNeuronsPerLayer; + auto numLayers_ = numLayers; + auto numHiddenLayers_ = numHiddenLayers; + auto AF_ = AF; + auto weights_ = weights; + auto bias_ = bias; + + auto calc_nn_op = KOKKOS_LAMBDA( const int atomindex ) + { + int attype = type( atomindex ); + // set input layer of NN + int layer_0, layer_lminusone; + layer_0 = (int)numSFperElem_( attype ); + + for ( int k = 0; k < layer_0; ++k ) + NN( atomindex, 0, k ) = G( atomindex, k ); + // forward propagation + for ( int l = 1; l < numLayers_; l++ ) + { + if ( l == 1 ) + layer_lminusone = layer_0; + else + layer_lminusone = numNeuronsPerLayer_( l - 1 ); + double dtmp; + for ( int i = 0; i < numNeuronsPerLayer_( l ); i++ ) + { + dtmp = 0.0; + for ( int j = 0; j < layer_lminusone; j++ ) + dtmp += weights_( attype, l - 1, i, j ) * + NN( atomindex, l - 1, j ); + dtmp += bias_( attype, l - 1, i ); + if ( AF_( l ) == 0 ) + { + NN( atomindex, l, i ) = dtmp; + dfdx( atomindex, l, i ) = 1.0; + } + else if ( AF_( l ) == 1 ) + { + dtmp = tanh( dtmp ); + NN( atomindex, l, i ) = dtmp; + dfdx( atomindex, l, i ) = 1.0 - dtmp * dtmp; + } + } + } + + E( atomindex ) = NN( atomindex, numLayers_ - 1, 0 ); + + // derivative of network w.r.t NN inputs + for ( int k = 0; k < numNeuronsPerLayer_( 0 ); k++ ) + { + for ( int i = 0; i < numNeuronsPerLayer_( 1 ); i++ ) + inner( atomindex, 0, i ) = + weights_( attype, 0, i, k ) * dfdx( atomindex, 1, i ); + + for ( int l = 1; l < numHiddenLayers_ + 1; l++ ) + { + for ( int i2 = 0; i2 < numNeuronsPerLayer_( l + 1 ); i2++ ) + { + outer( atomindex, l - 1, i2 ) = 0.0; + + for ( int i1 = 0; i1 < numNeuronsPerLayer_( l ); i1++ ) + outer( atomindex, l - 1, i2 ) += + weights_( attype, l, i2, i1 ) * + inner( atomindex, l - 1, i1 ); + outer( atomindex, l - 1, i2 ) *= + dfdx( atomindex, l + 1, i2 ); + + if ( l < numHiddenLayers_ ) + inner( atomindex, l, i2 ) = + outer( atomindex, l - 1, i2 ); + } + } + dEdG( atomindex, k ) = outer( atomindex, numHiddenLayers_ - 1, 0 ); + } + }; + Kokkos::parallel_for( "Mode::calculateAtomicNeuralNetworks", policy, + calc_nn_op ); + Kokkos::fence(); +} + +template +template +void ModeKokkos::calculateForces( + 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; + + 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 ) + { + double pfcij = 0.0; + double pdfcij = 0.0; + double rij, r2ij; + T_F_FLOAT dxij, dyij, dzij; + 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 ); + ++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_ ); + + for(int jj = 0; jj < num_neighs_i; jj++) + { + 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 ) + { + // 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( "Mode::calculateRadialForces", + policy, calc_radial_force_op ); + Kokkos::fence(); + + auto calc_angular_force_op = KOKKOS_LAMBDA( const int i ) + { + 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; + + 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 ) + { + 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_ ); + + for(int jj = 0; jj < num_neighs_i; jj++) + { + 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 ) + { + // 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++) + { + T_INT k = neighs_j(kk); + 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. + 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( "Mode::calculateAngularForces", + policy, calc_angular_force_op ); + Kokkos::fence(); + + return; +} + +} +#endif diff --git a/src/libnnpif/CabanaMD/typesCabana.h b/src/libnnpif/Kokkos/typesKokkos.h similarity index 95% rename from src/libnnpif/CabanaMD/typesCabana.h rename to src/libnnpif/Kokkos/typesKokkos.h index 66d32bd2d..a177ef023 100644 --- a/src/libnnpif/CabanaMD/typesCabana.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; 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 ###############################################################################