ver 2021.9.28

USER-PACE/ace_abstract_basis.cpp → ML-PACE/ace_abstract_basis.cpp

@@ -28,6 +28,7 @@
 // Created by Lysogorskiy Yury on 28.04.2020.
 
 #include "ace_abstract_basis.h"
+#include "ace_radial.h"
 
 ////embedding function
 ////case nemb = 1 only implementation
@@ -93,21 +94,27 @@ void ACEAbstractBasisSet::inner_cutoff(DOUBLE_TYPE rho_core, DOUBLE_TYPE rho_cut
         fcut = 1;
         dfcut = 0;
     } else {
-        fcut = 0.5 * (1 + cos(M_PI * (rho_core - rho_low) / drho_cut));
-        dfcut = -0.5 * sin(M_PI * (rho_core - rho_low) / drho_cut) * M_PI / drho_cut;
+        cutoff_func_poly(rho_core, rho_cut, drho_cut, fcut, dfcut);
     }
 }
 
 void ACEAbstractBasisSet::FS_values_and_derivatives(Array1D<DOUBLE_TYPE> &rhos, DOUBLE_TYPE &value,
-                                                    Array1D<DOUBLE_TYPE> &derivatives, DENSITY_TYPE ndensity) {
+                                                    Array1D<DOUBLE_TYPE> &derivatives,
+                                                    SPECIES_TYPE mu_i
+) {
     DOUBLE_TYPE F, DF = 0, wpre, mexp;
+    DENSITY_TYPE ndensity = map_embedding_specifications.at(mu_i).ndensity;
     for (int p = 0; p < ndensity; p++) {
-        wpre = FS_parameters.at(p * ndensity + 0);
-        mexp = FS_parameters.at(p * ndensity + 1);
-        if (this->npoti == "FinnisSinclair")
+
+        wpre = map_embedding_specifications.at(mu_i).FS_parameters.at(p * 2 + 0);
+        mexp = map_embedding_specifications.at(mu_i).FS_parameters.at(p * 2 + 1);
+        string npoti = map_embedding_specifications.at(mu_i).npoti;
+
+        if (npoti == "FinnisSinclair")
             Fexp(rhos(p), mexp, F, DF);
-        else if (this->npoti == "FinnisSinclairShiftedScaled")
+        else if (npoti == "FinnisSinclairShiftedScaled")
             FexpShiftedScaled(rhos(p), mexp, F, DF);
+
         value += F * wpre; // * weight (wpre)
         derivatives(p) = DF * wpre;// * weight (wpre)
     }
@@ -144,8 +151,9 @@ ACEAbstractBasisSet::~ACEAbstractBasisSet() {
 
 void ACEAbstractBasisSet::_copy_scalar_memory(const ACEAbstractBasisSet &src) {
     deltaSplineBins = src.deltaSplineBins;
-    FS_parameters = src.FS_parameters;
-    npoti = src.npoti;
+
+    map_bond_specifications = src.map_bond_specifications;
+    map_embedding_specifications = src.map_embedding_specifications;
 
     nelements = src.nelements;
     rankmax = src.rankmax;
@@ -157,10 +165,6 @@ void ACEAbstractBasisSet::_copy_scalar_memory(const ACEAbstractBasisSet &src) {
 
     spherical_harmonics = src.spherical_harmonics;
 
-    rho_core_cutoffs = src.rho_core_cutoffs;
-    drho_core_cutoffs = src.drho_core_cutoffs;
-
-
     E0vals = src.E0vals;
 }
 
@@ -183,3 +187,50 @@ SPECIES_TYPE ACEAbstractBasisSet::get_species_index_by_name(const string &elemna
     return -1;
 }
 
+bool compare(const ACEAbstractBasisFunction &b1, const ACEAbstractBasisFunction &b2) {
+    if (b1.rank < b2.rank) return true;
+    else if (b1.rank > b2.rank) return false;
+    //now b1.rank==b2.rank
+
+
+//    if (b1.rankL < b2.rankL) return true;
+//    else if (b1.rankL > b2.rankL) return false;
+    //now b1.rankL==b2.rankL
+
+    if (b1.mu0 < b2.mu0) return true;
+    else if (b1.mu0 > b2.mu0) return false;
+
+    for (RANK_TYPE r = 0; r < b1.rank; r++)
+        if (b1.mus[r] < b2.mus[r]) return true;
+        else if (b1.mus[r] > b2.mus[r]) return false;
+    // now mus are checked to be equal
+
+    for (RANK_TYPE r = 0; r < b1.rank; r++)
+        if (b1.ns[r] < b2.ns[r]) return true;
+        else if (b1.ns[r] > b2.ns[r]) return false;
+    // now ns are checked to be equal
+
+    for (RANK_TYPE r = 0; r < b1.rank; r++)
+        if (b1.ls[r] < b2.ls[r]) return true;
+        else if (b1.ls[r] > b2.ls[r]) return false;
+    // now ls are checked to be equal
+
+//    for (RANK_TYPE r = 0; r < b1.rankL; r++)
+//        if (b1.LS[r] < b2.LS[r]) return true;
+//        else if (b1.LS[r] > b2.LS[r]) return false;
+    // now LS are checked to be equal
+
+    if (b1.num_ms_combs < b2.num_ms_combs) return true;
+    else if (b1.num_ms_combs > b2.num_ms_combs) return false;
+    // now num_ms_combs are checked to be equal
+
+    // ms_combs:
+    // size =  num_ms_combs * rank,
+    // effective shape: [num_ms_combs][rank]
+    for (SHORT_INT_TYPE i = 0; i < b1.num_ms_combs * b1.rank; i++)
+        if (b1.ms_combs[i] < b2.ms_combs[i]) return true;
+        else if (b1.ms_combs[i] > b2.ms_combs[i]) return false;
+
+    return false; // ! b1<b2
+}
+

USER-PACE/ace_abstract_basis.h → ML-PACE/ace_abstract_basis.h

@@ -33,6 +33,8 @@
 
 #include <vector>
 #include <string>
+#include <map>
+#include <tuple>
 
 #include "ace_c_basisfunction.h"
 #include "ace_contigous_array.h"
@@ -42,6 +44,138 @@
 
 using namespace std;
 
+struct ACEEmbeddingSpecification {
+    DENSITY_TYPE ndensity;
+    vector<DOUBLE_TYPE> FS_parameters; ///< parameters for cluster functional, see Eq.(3) in implementation notes or Eq.(53) in <A HREF="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.014104">  PRB 99, 014104 (2019) </A>
+    string npoti = "FinnisSinclair"; ///< FS and embedding function combination
+
+    DOUBLE_TYPE rho_core_cutoff;
+    DOUBLE_TYPE drho_core_cutoff;
+
+    //TODO: update method
+    string to_string() const {
+        stringstream ss;
+        ss << "ACEEmbeddingSpecification(npoti=" << npoti << ", ndensity=" << ndensity << ", ";
+        ss << "FS_parameter=(";
+        for (const auto &p: FS_parameters)
+            ss << p << ",";
+        ss << "))";
+        return ss.str();
+    }
+
+    YAML::Node to_YAML() const {
+        YAML::Node emb_yaml;
+        emb_yaml.SetStyle(YAML::EmitterStyle::Flow);
+        emb_yaml["ndensity"] = (int) this->ndensity;
+        emb_yaml["FS_parameters"] = this->FS_parameters;
+
+        emb_yaml["npoti"] = this->npoti;
+        emb_yaml["rho_core_cutoff"] = this->rho_core_cutoff;
+        emb_yaml["drho_core_cutoff"] = this->drho_core_cutoff;
+
+        return emb_yaml;
+    }
+};
+
+struct ACEBondSpecification {
+    NS_TYPE nradmax;
+    LS_TYPE lmax;
+    NS_TYPE nradbasemax;
+
+    string radbasename;
+
+    std::vector<DOUBLE_TYPE> radparameters; // i.e. lambda
+    vector<vector<vector<DOUBLE_TYPE>>> radcoefficients;///< crad_nlk order: [n=0..nradmax-1][l=0..lmax][k=0..nradbase-1]
+
+    //hard-core repulsion
+    DOUBLE_TYPE prehc;
+    DOUBLE_TYPE lambdahc;
+
+    DOUBLE_TYPE rcut;
+    DOUBLE_TYPE dcut;
+
+    //inner cutoff
+    DOUBLE_TYPE rcut_in = 0;
+    DOUBLE_TYPE dcut_in = 0;
+    string inner_cutoff_type = "distance"; // new behaviour is default
+
+    bool operator==(const ACEBondSpecification &another) const {
+        return (nradbasemax == another.nradbasemax) && (lmax == another.lmax) &&
+               (nradbasemax == another.nradbasemax) && (radbasename == another.radbasename) &&
+               (radparameters == another.radparameters) && (radcoefficients == another.radcoefficients) &&
+               (prehc == another.prehc) && (lambdahc == another.lambdahc) && (rcut == another.rcut) &&
+               (dcut == another.dcut) && (rcut_in == another.rcut_in) && (dcut_in == another.dcut_in) &&
+               (inner_cutoff_type == another.inner_cutoff_type);
+    }
+
+    bool operator!=(const ACEBondSpecification &another) const {
+        return !(*this == (another));
+    }
+
+    string to_string() const {
+        stringstream ss;
+        ss << "ACEBondSpecification(nradmax=" << nradmax << ", lmax=" << lmax << ", nradbasemax=" << nradbasemax
+           << ", radbasename=" <<
+           radbasename << ", crad=(" << radcoefficients.at(0).at(0).at(0) << "...), ";
+        ss << "rcut=" << rcut << ", dcut=" << dcut;
+        ss << ", rcut_in=" << rcut_in << ", dcut_in=" << dcut_in;
+        ss << ", inner_cutoff_type=" << inner_cutoff_type;
+        if (prehc > 0)
+            ss << ", core-rep: [prehc=" << prehc << ", lambdahc=" << lambdahc << "]";
+        ss << ")";
+        return ss.str();
+    }
+
+    void from_YAML(YAML::Node bond_yaml) {
+        radbasename = bond_yaml["radbasename"].as<string>();
+
+        if(radbasename=="ACE.jl.base") {
+
+        } else {
+            nradmax = bond_yaml["nradmax"].as<NS_TYPE>();
+            lmax = bond_yaml["lmax"].as<LS_TYPE>();
+            nradbasemax = bond_yaml["nradbasemax"].as<NS_TYPE>();
+            radparameters = bond_yaml["radparameters"].as<vector<DOUBLE_TYPE>>();
+            radcoefficients = bond_yaml["radcoefficients"].as<vector<vector<vector<DOUBLE_TYPE>>>>();
+            prehc = bond_yaml["prehc"].as<DOUBLE_TYPE>();
+            lambdahc = bond_yaml["lambdahc"].as<DOUBLE_TYPE>();
+            rcut = bond_yaml["rcut"].as<DOUBLE_TYPE>();
+            dcut = bond_yaml["dcut"].as<DOUBLE_TYPE>();
+
+            if (bond_yaml["rcut_in"]) rcut_in = bond_yaml["rcut_in"].as<DOUBLE_TYPE>();
+            if (bond_yaml["dcut_in"]) dcut_in = bond_yaml["dcut_in"].as<DOUBLE_TYPE>();
+            if (bond_yaml["inner_cutoff_type"])
+                inner_cutoff_type = bond_yaml["inner_cutoff_type"].as<string>();
+            else
+                inner_cutoff_type = "density"; // default value to read for backward compatibility
+        }
+    }
+
+    YAML::Node to_YAML() const {
+        YAML::Node bond_yaml;
+        bond_yaml.SetStyle(YAML::EmitterStyle::Flow);
+        bond_yaml["nradmax"] = (int) this->nradmax;
+        bond_yaml["lmax"] = (int) this->lmax;
+        bond_yaml["nradbasemax"] = (int) this->nradbasemax;
+
+        bond_yaml["radbasename"] = this->radbasename;
+        bond_yaml["radparameters"] = this->radparameters;
+        bond_yaml["radcoefficients"] = this->radcoefficients;
+
+        bond_yaml["prehc"] = this->prehc;
+        bond_yaml["lambdahc"] = this->lambdahc;
+        bond_yaml["rcut"] = this->rcut;
+        bond_yaml["dcut"] = this->dcut;
+
+        bond_yaml["rcut_in"] = this->rcut_in;
+        bond_yaml["dcut_in"] = this->dcut_in;
+        bond_yaml["inner_cutoff_type"] = this->inner_cutoff_type;
+
+        return bond_yaml;
+    }
+};
+
+
 /**
  * Abstract basis set class
  */
@@ -56,18 +190,18 @@ class ACEAbstractBasisSet {
     DOUBLE_TYPE cutoffmax = 0;  ///< maximum value of cutoff distance among all species in basis set
     DOUBLE_TYPE deltaSplineBins = 0;  ///< Spline interpolation density
 
-    string npoti = "FinnisSinclair"; ///< FS and embedding function combination
-
     string *elements_name = nullptr; ///< Array of elements name for mapping from index (0..nelements-1) to element symbol (string)
+    map<string, SPECIES_TYPE> elements_to_index_map;
+
 
     AbstractRadialBasis *radial_functions = nullptr; ///< object to work with radial functions
     ACECartesianSphericalHarmonics spherical_harmonics; ///< object to work with spherical harmonics in Cartesian representation
 
+    bool is_sort_functions = true; ///< flag to specify, whether to sort the basis functions or preserve the order
+    //for multispecies
 
-    Array1D<DOUBLE_TYPE> rho_core_cutoffs; ///< energy-based inner cut-off
-    Array1D<DOUBLE_TYPE> drho_core_cutoffs; ///< decay of energy-based inner cut-off
-
-    vector<DOUBLE_TYPE> FS_parameters;  ///< parameters for cluster functional, see Eq.(3) in implementation notes or Eq.(53) in <A HREF="https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.014104">  PRB 99, 014104 (2019) </A>
+    map<SPECIES_TYPE, ACEEmbeddingSpecification> map_embedding_specifications;
+    map<pair<SPECIES_TYPE, SPECIES_TYPE>, ACEBondSpecification> map_bond_specifications;
 
     // E0 values
     Array1D<DOUBLE_TYPE> E0vals;
@@ -106,7 +240,7 @@ class ACEAbstractBasisSet {
      * @param ndensity  number \f$ P \f$ of densities to use
      */
     void FS_values_and_derivatives(Array1D<DOUBLE_TYPE> &rhos, DOUBLE_TYPE &value, Array1D<DOUBLE_TYPE> &derivatives,
-                                   DENSITY_TYPE ndensity);
+                                   SPECIES_TYPE mu_i);
 
     /**
      * Computing hard core pairwise repulsive potential \f$ f_{cut}(\rho_i^{(\textrm{core})})\f$ and its derivative,
@@ -160,11 +294,19 @@ class ACEAbstractBasisSet {
      * @param src source object to copy from
      */
     virtual void _copy_scalar_memory(const ACEAbstractBasisSet &src);
+
+    virtual vector<DOUBLE_TYPE> get_all_coeffs() const = 0;
+
+    virtual vector<vector<SPECIES_TYPE>> get_all_coeffs_mask() const = 0;
+
+    virtual void set_all_coeffs(const vector<DOUBLE_TYPE> &coeffs) = 0;
 };
 
 void Fexp(DOUBLE_TYPE rho, DOUBLE_TYPE mexp, DOUBLE_TYPE &F, DOUBLE_TYPE &DF);
 
 void FexpShiftedScaled(DOUBLE_TYPE rho, DOUBLE_TYPE mexp, DOUBLE_TYPE &F, DOUBLE_TYPE &DF);
 
+bool compare(const ACEAbstractBasisFunction &b1, const ACEAbstractBasisFunction &b2);
+
 #endif //ACE_EVALUATOR_ACE_ABSTRACT_BASIS_H
 

USER-PACE/ace_array2dlm.h → ML-PACE/ace_array2dlm.h

@@ -492,6 +492,7 @@ class Array4DLM : public ContiguousArrayND<T> {
     }
 
 #ifdef MULTIARRAY_INDICES_CHECK
+
     /**
      * Check if indices (i0, l,m) are within array
      */
@@ -519,10 +520,12 @@ class Array4DLM : public ContiguousArrayND<T> {
 
         size_t ii = i0 * s[0] + i1 * s[1] + l * (l + 1) + m;
         if (ii >= size) {
-            fprintf(stderr, "%s: index = %ld out of range %ld\n", array_name.c_str(), ii, size);
+            fprintf(stderr, "%s: index = %ld (i0=%ld, i1=%ld, l=%d, m=%d) out of range %ld\n", array_name.c_str(), ii, i0,
+                    i1, l, m, size);
             exit(EXIT_FAILURE);
         }
     }
+
 #endif
 
     /**