change ships_radial to acejl_radial with splines.

ML-PACE/ace-evaluator/ace_c_basis.cpp

@@ -29,7 +29,7 @@
 #include <fstream>
 
 #include "ace-evaluator/ace_c_basis.h"
-#include "ace-evaluator/ships_radial.h"
+#include "ace-evaluator/acejl_radial.h"
 
 using namespace std;
 
@@ -679,7 +679,7 @@ void ACECTildeBasisSet::load(const string filename) {
     if ((radbasename == "ChebExpCos") || (radbasename == "ChebPow")) {
         _load_radial_ACERadial(fptr, filename, radbasename);
     } else if (radbasename == "ACE.jl.Basic") {
-        _load_radial_SHIPsBasic(fptr, filename, radbasename);
+        _load_radial_ACEjlBasic(fptr, filename, radbasename);
     } else {
         throw invalid_argument(
                 ("File '" + filename + "': I don't know how to read radbasename = " + radbasename).c_str());
@@ -1018,54 +1018,10 @@ void ACECTildeBasisSet::_load_radial_ACERadial(FILE *fptr,
         }
 }
 
-void ACECTildeBasisSet::_load_radial_SHIPsBasic(FILE *fptr,
+void ACECTildeBasisSet::_load_radial_ACEjlBasic(FILE *fptr,
                                                 const string filename,
                                                 const string radbasename) {
-    // create a radial basis object, and read it from the file pointer
-    SHIPsRadialFunctions *ships_radial_functions = new SHIPsRadialFunctions();
-
-    ships_radial_functions->nelements = nelements;
-    ships_radial_functions->radbasis.init(nelements, nelements, "SHIPsRadialFunctions::radbasis");
-    ships_radial_functions->fread(fptr);
-
-    _post_load_radial_SHIPsBasic(ships_radial_functions);
-}
-
-void ACECTildeBasisSet::_post_load_radial_SHIPsBasic(
-        SHIPsRadialFunctions *ships_radial_functions) {//mimic ships_radial_functions to ACERadialFunctions
-    ships_radial_functions->nradial = ships_radial_functions->get_maxn();
-    ships_radial_functions->nradbase = ships_radial_functions->get_maxn();
-
-    nradbase = ships_radial_functions->get_maxn();
-    nradmax = ships_radial_functions->get_maxn();
-    cutoffmax = ships_radial_functions->get_rcut();
-    deltaSplineBins = 0.001;
-
-    ships_radial_functions->nradbase = nradbase;
-    ships_radial_functions->lmax = lmax;
-    ships_radial_functions->nradial = nradmax;
-    ships_radial_functions->gr.init(nradbase, "gr");
-    ships_radial_functions->dgr.init(nradbase, "dgr");
-
-    ships_radial_functions->fr.init(nradmax, lmax + 1, "fr");
-    ships_radial_functions->dfr.init(nradmax, lmax + 1, "dfr");
-    ships_radial_functions->crad.init(nelements, nelements, nradmax, (lmax + 1), nradbase, "crad");
-    ships_radial_functions->crad.fill(0.);
-
-
-    if (radial_functions) delete radial_functions;
-    radial_functions = ships_radial_functions;
-    radial_functions->prehc.fill(0);
-    radial_functions->lambdahc.fill(1);
-    radial_functions->lambda.fill(0);
-
-    radial_functions->cut.init(nelements, nelements, "cut");
-    radial_functions->dcut.init(nelements, nelements, "dcut");
-
-    radial_functions->cut.fill(ships_radial_functions->get_rcut());
-    radial_functions->dcut.fill(0);
-
-    radial_functions->crad.fill(0);
+    throw runtime_error("_load_radial_ACEjlBasic is deprecated");
 }
 
 vector<vector<SPECIES_TYPE>> ACECTildeBasisSet::get_all_coeffs_mask() const {
@@ -1244,18 +1200,18 @@ void ACECTildeBasisSet::load_yaml(const string &yaml_file_name) {
     this->cutoffmax = 0;
 
     // check, if bonds::[]::radbasename=="ACE.jl.radbase"
-    bool ACE_jl_radbase = false;
+    bool ACEjl_radbase = false;
     bool PACE_radbase = false;
     for (const auto &p: yaml_map_bond_specifications) {
         auto bond_yaml = p.second;
         string radbasename = bond_yaml["radbasename"].as<string>();
         if (radbasename.rfind("ACE.jl", 0) == 0)
-            ACE_jl_radbase = true;
+            ACEjl_radbase = true;
         else
             PACE_radbase = true;
     }
     // check if both type of radbase -> inconsistency
-    if (ACE_jl_radbase & PACE_radbase) {
+    if (ACEjl_radbase & PACE_radbase) {
         throw invalid_argument(
                 "Only ACE.jl.* or PACE's radial basis are possible, but both types are used simultaneously.");
     }
@@ -1326,19 +1282,41 @@ void ACECTildeBasisSet::load_yaml(const string &yaml_file_name) {
             }
         }
         ///////////////////////////////////////////////////////////////////
-    } else if (ACE_jl_radbase) {
+    } else if (ACEjl_radbase) {
         ///////////////////////////////////////////////////////////////////
-        //read  lmax from YACE
-        if (ctilde_basis_yaml["lmax"])
-            this->lmax = ctilde_basis_yaml["lmax"].as<LS_TYPE>();
-        else
-            throw invalid_argument(
-                    "For `ACE.jl.*` radbase functions, `lmax` should be provided in the YACE separately.");
-        // no need to store map_bond_specifications, only SHIPsRadialFunctions
-        SHIPsRadialFunctions *ships_radial_functions = new SHIPsRadialFunctions();
-        ships_radial_functions->init(nelements);
-        ships_radial_functions->read_yaml(ctilde_basis_yaml);
-        _post_load_radial_SHIPsBasic(ships_radial_functions);
+
+        lmax = ctilde_basis_yaml["lmax"].as<LS_TYPE>();
+
+        // loop through bonds once to get what is needed for `init`
+        // (remaining information is extracted later in read_yaml)
+        vector<vector<string>> radbasename_ij(nelements, vector<string>(nelements));
+        cutoffmax = 0.0;
+        nradmax = 0;
+        nradbase = 0;
+        for (const auto &p: yaml_map_bond_specifications) {
+            pair<SPECIES_TYPE, SPECIES_TYPE> bond_pair = make_pair(p.first[0], p.first[1]);
+            auto bond_yaml = p.second;
+            radbasename_ij.at(bond_pair.first).at(bond_pair.second) = bond_yaml["radbasename"].as<string>();
+            if (bond_yaml["rcut"].as<DOUBLE_TYPE>() > cutoffmax)
+                cutoffmax = bond_yaml["rcut"].as<DOUBLE_TYPE>();
+            if (bond_yaml["nradial"].as<NS_TYPE>() > nradmax) {
+                nradmax = bond_yaml["nradial"].as<NS_TYPE>();
+                nradbase = bond_yaml["nradial"].as<NS_TYPE>();
+            }
+        }
+
+        // create and initialize acejl_radial_functions
+        ACEjlRadialFunctions* acejl_radial_functions = new ACEjlRadialFunctions();
+        acejl_radial_functions->init(nradbase, lmax, nradmax,
+                                     -1.0, nelements, radbasename_ij);
+
+        // read remaining information from yaml
+        acejl_radial_functions->read_yaml(ctilde_basis_yaml);
+
+        // make radial_functions point to the acejl basis
+        if (radial_functions) delete radial_functions;
+        radial_functions = acejl_radial_functions;
+
     }
     ///////////////////////////////////////////////////////////////////
 

ML-PACE/ace-evaluator/ace_c_basis.h

@@ -31,7 +31,7 @@
 #define ACE_C_BASIS_H
 
 #include "ace-evaluator/ace_flatten_basis.h"
-#include "ace-evaluator/ships_radial.h"
+#include "ace-evaluator/acejl_radial.h"
 
 typedef vector<vector<ACECTildeBasisFunction>> C_tilde_full_basis_vector2d;
 
@@ -129,7 +129,7 @@ class ACECTildeBasisSet : public ACEFlattenBasisSet {
                                 const string filename,
                                 const string radbasename);
 
-    void _load_radial_SHIPsBasic(FILE *fptr,
+    void _load_radial_ACEjlBasic(FILE *fptr,
                                  const string filename,
                                  const string radbasename);
 
@@ -168,8 +168,6 @@ class ACECTildeBasisSet : public ACEFlattenBasisSet {
 
     void set_all_coeffs(const vector<DOUBLE_TYPE> &coeffs) override;
 
-
-    void _post_load_radial_SHIPsBasic(SHIPsRadialFunctions *ships_radial_functions);
 };
 
 #endif //ACE_C_BASIS_H

ML-PACE/ace-evaluator/acejl_radial.cpp

@@ -0,0 +1,172 @@
+/*
+ * Performant implementation of atomic cluster expansion and interface to LAMMPS
+ *
+ * Copyright 2021  (c) Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1,
+ * Sarath Menon^1, Matteo Rinaldi^1, Thomas Hammerschmidt^1, Matous Mrovec^1,
+ * Aidan Thompson^3, Gabor Csanyi^2, Christoph Ortner^4, Ralf Drautz^1
+ *
+ * ^1: Ruhr-University Bochum, Bochum, Germany
+ * ^2: University of Cambridge, Cambridge, United Kingdom
+ * ^3: Sandia National Laboratories, Albuquerque, New Mexico, USA
+ * ^4: University of British Columbia, Vancouver, BC, Canada
+ *
+ *
+ * See the LICENSE file.
+ * This FILENAME 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 <http://www.gnu.org/licenses/>.
+ */
+
+// Created by Christoph Ortner  on 03.06.2020
+// Adapted to use splines by Chuck Witt (April 2023)
+
+#include "ace-evaluator/acejl_radial.h"
+
+#include <functional>
+#include <cmath>
+#include <string>
+
+using namespace std;
+
+void ACEjlRadialFunctions::read_yaml(YAML_PACE::Node node) {
+
+    auto yaml_map_bond_specifications = node["bonds"].as<map<vector<int>, YAML_PACE::Node>>();
+
+    // loop over element pairs
+    for (const auto &p: yaml_map_bond_specifications) {
+
+        SPECIES_TYPE mu_i = p.first[0];
+        SPECIES_TYPE mu_j = p.first[1];
+
+        // extract spline information for this pair 
+        auto bond_yaml = p.second;
+        cut(mu_i,mu_j) = bond_yaml["rcut"].as<DOUBLE_TYPE>();
+        int nbins = bond_yaml["nbins"].as<int>();
+        auto splinenodalvals = bond_yaml["splinenodalvals"].as<map<int,vector<DOUBLE_TYPE>>>();
+        auto splinenodalderivs = bond_yaml["splinenodalderivs"].as<map<int,vector<DOUBLE_TYPE>>>();
+
+        // set up spline interpolator (see SplineInterpolator::setupSplines)
+        // note that for SplineInterpolators, ntot=nlut is defined such 
+        // that if r*rscalelookup>=ntot then f(r/delta)=0
+        SplineInterpolator& spl = splines(mu_i, mu_j);
+        spl.cutoff = cut(mu_i,mu_j);
+        spl.ntot = nbins;
+        spl.deltaSplineBins = spl.cutoff / spl.ntot;
+        spl.num_of_functions = splinenodalvals.size();
+        spl.values.resize(spl.num_of_functions);
+        spl.derivatives.resize(spl.num_of_functions);
+        spl.second_derivatives.resize(spl.num_of_functions);
+        spl.nlut = spl.ntot;
+        spl.rscalelookup = (DOUBLE_TYPE) spl.nlut / spl.cutoff;
+        spl.invrscalelookup = 1.0 / spl.rscalelookup;
+        spl.lookupTable.init(spl.ntot+1, spl.num_of_functions, 4);
+        for (int n=0; n<spl.nlut; ++n) {
+            for (int func_id=0; func_id<spl.num_of_functions; func_id++) {
+                const auto d = spl.deltaSplineBins;
+                DOUBLE_TYPE f0 = splinenodalvals[func_id][n];
+                DOUBLE_TYPE f1 = splinenodalvals[func_id][n+1];
+                DOUBLE_TYPE f0d1 = splinenodalderivs[func_id][n]*d;
+                DOUBLE_TYPE f1d1 = splinenodalderivs[func_id][n+1]*d;
+                // store c0 + c1*r + c2*r^2 + c3*r^3 coefficients
+                spl.lookupTable(n,func_id,0) = f0;
+                spl.lookupTable(n,func_id,1) = f0d1;
+                spl.lookupTable(n,func_id,2) = 3.0*(f1-f0)-f1d1-2.0*f0d1;
+                spl.lookupTable(n,func_id,3) = -2.0*(f1-f0)+f1d1+f0d1;
+            }
+        }
+        // haven't yet set the last element in the lookup table; set to zero
+        for (int func_id=0; func_id<spl.num_of_functions; func_id++) {
+            spl.lookupTable(spl.nlut,func_id,0) = 0.0;
+            spl.lookupTable(spl.nlut,func_id,1) = 0.0;
+            spl.lookupTable(spl.nlut,func_id,2) = 0.0;
+            spl.lookupTable(spl.nlut,func_id,3) = 0.0;
+        }
+    }
+}
+
+void ACEjlRadialFunctions::setuplookupRadspline() {
+}
+
+void ACEjlRadialFunctions::init(
+        NS_TYPE nradb, LS_TYPE lmax, NS_TYPE nradial,
+        DOUBLE_TYPE deltaSplineBins, SPECIES_TYPE nelements,
+        vector<vector<string>> radbasename) {
+
+    // these spline members are unique to ACEjlRadialFunctions
+    splines.init(nelements, nelements, "splines");
+
+    // ... and the rest are inherited from AbstractRadialBasis
+    this->nelements = nelements;
+    cut.init(nelements, nelements, "cut");
+    cut.fill(0.0);  // zeroed here, updated later
+    dcut.init(nelements, nelements, "dcut");
+    dcut.fill(1.0);
+
+    this->deltaSplineBins = deltaSplineBins;
+    this->lmax = lmax;
+    this->nradial = nradial;
+    nradbase = nradb;
+
+    radbasenameij = radbasename;
+
+    gr.init(nradb, "gr");
+    dgr.init(nradb, "dgr");
+    d2gr.init(nradbase, "d2gr");
+
+    fr.init(nradial, lmax + 1, "fr");
+    dfr.init(nradial, lmax + 1, "dfr");
+    d2fr.init(nradial, lmax + 1, "d2fr");
+
+    cr = 0.0;
+    dcr = 0.0;
+
+    crad.init(nelements, nelements, nradial, (lmax + 1), nradb, "crad");
+    crad.fill(0.0);
+
+    cut_in.init(nelements, nelements, "cut_in");
+    cut_in.fill(0.0);
+
+    dcut_in.init(nelements, nelements, "dcut_in");
+    dcut_in.fill(1e-5);
+
+    lambda.init(nelements, nelements, "lambda");
+    lambda.fill(1.0);
+
+    prehc.init(nelements, nelements, "prehc");
+    prehc.fill(0.0);
+    lambdahc.init(nelements, nelements, "lambdahc");
+    lambdahc.fill(1.0);
+}
+
+void ACEjlRadialFunctions::evaluate(
+        DOUBLE_TYPE r, NS_TYPE nradbase_c, NS_TYPE nradial_c,
+        SPECIES_TYPE mu_i, SPECIES_TYPE mu_j,
+        bool calc_second_derivatives) {
+
+    auto &spline = splines(mu_i, mu_j);
+    spline.calcSplines(r, calc_second_derivatives);
+
+    for (NS_TYPE n = 0; n < gr.get_dim(0); n++) {
+        gr(n) = spline.values(n);
+        dgr(n) = spline.derivatives(n);
+        if (calc_second_derivatives)
+            d2gr(n) = spline.second_derivatives(n);
+    }
+
+    for (NS_TYPE n = 0; n < fr.get_dim(0); n++) {
+        for (LS_TYPE l = 0; l < fr.get_dim(1); l++) {
+            fr(n,l) = spline.values(n);
+            dfr(n,l) = spline.derivatives(n);
+            if (calc_second_derivatives)
+                d2fr(n,l) = spline.second_derivatives(n);
+        }
+    }
+}