Making MLIP as a submodule

CMakeLists.txt

+cmake_minimum_required(VERSION 3.12)
+
+project(Tadah.MLIP)
+
+file(GLOB_RECURSE MLIP_SRC2 CONFIGURE_DEPENDS "*.cpp")
+set(MLIP_SRC ${MLIP_SRC2} PARENT_SCOPE)
+
+
+#add_library(tadah.mlip
+#    ${MLIP_SRC}
+#    )
+

analytics/analytics.cpp

+#include "analytics.h"
+#include "statistics.h"
+
+Analytics::Analytics(const StructureDB &st, const StructureDB &stp):
+    st(st),
+    stp(stp)
+{
+    if (st.size()!=stp.size())
+        throw std::runtime_error("Containers differ in size.");
+}
+
+t_type Analytics::calc_e_mae() const{
+
+    t_type emae_vec(st.dbidx.size()-1);
+    double emae=0;
+    size_t dbidx=0;
+    size_t N=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        // assume stp(i).natoms()==st(i).natoms()
+        emae += std::abs(st(i).energy - stp(i).energy)/st(i).natoms();
+        N++;
+        if (i+1==st.dbidx[dbidx+1]) {
+            emae_vec(dbidx)=emae/N;
+            emae=0;
+            N=0;
+            dbidx++;
+        }
+    }
+    return emae_vec;
+}
+
+t_type Analytics::calc_f_mae() const{
+
+    t_type fmae_vec(st.dbidx.size()-1);
+    double fmae=0;
+    size_t dbidx=0;
+    size_t N=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        // assume stp(i).natoms()==st(i).natoms()
+
+        for (size_t a=0; a<st(i).natoms(); ++a) {
+            for (size_t k=0; k<3; ++k) {
+                fmae += std::abs(st(i).atoms[a].force(k) - stp(i).atoms[a].force(k));
+                N++;
+            }
+        }
+
+        if (i+1==st.dbidx[dbidx+1]) {
+            fmae_vec(dbidx)=fmae/N;
+            fmae=0;
+            N=0;
+            dbidx++;
+        }
+    }
+    return fmae_vec;
+}
+
+t_type Analytics::calc_s_mae() const {
+
+    t_type smae_vec(st.dbidx.size()-1);
+    double smae=0;
+    size_t dbidx=0;
+    size_t N=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        for (size_t x=0; x<3; ++x) {
+            for (size_t y=x; y<3; ++y) {
+                smae += std::abs(st(i).stress(x,y) - stp(i).stress(x,y));
+                N++;
+            }
+        }
+        if (i+1==st.dbidx[dbidx+1]) {
+            smae_vec(dbidx)=smae/N;
+            smae=0;
+            N=0;
+            dbidx++;
+        }
+    }
+    return smae_vec;
+}
+
+t_type Analytics::calc_e_rmse() const{
+
+    t_type ermse_vec(st.dbidx.size()-1);
+    double ermse=0;
+    size_t dbidx=0;
+    size_t N=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        // assume stp(i).natoms()==st(i).natoms()
+        ermse += std::pow((st(i).energy - stp(i).energy)/st(i).natoms(),2);
+        N++;
+        if (i+1==st.dbidx[dbidx+1]) {
+            ermse_vec(dbidx)=std::sqrt(ermse/N);
+            ermse=0;
+            N=0;
+            dbidx++;
+        }
+    }
+    return ermse_vec;
+}
+
+t_type Analytics::calc_f_rmse() const{
+
+    t_type frmse_vec(st.dbidx.size()-1);
+    double frmse=0;
+    size_t dbidx=0;
+    size_t N=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        // assume stp(i).natoms()==st(i).natoms()
+
+        for (size_t a=0; a<st(i).natoms(); ++a) {
+            for (size_t k=0; k<3; ++k) {
+                frmse += std::pow(st(i).atoms[a].force(k) - stp(i).atoms[a].force(k),2);
+                N++;
+            }
+        }
+
+        if (i+1==st.dbidx[dbidx+1]) {
+            frmse_vec(dbidx)=sqrt(frmse/N);
+            frmse=0;
+            N=0;
+            dbidx++;
+        }
+    }
+    return frmse_vec;
+}
+
+t_type Analytics::calc_s_rmse() const {
+
+    t_type srmse_vec(st.dbidx.size()-1);
+    double srmse=0;
+    size_t dbidx=0;
+    size_t N=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        for (size_t x=0; x<3; ++x) {
+            for (size_t y=x; y<3; ++y) {
+                srmse += std::pow(st(i).stress(x,y) - stp(i).stress(x,y),2);
+                N++;
+            }
+        }
+        if (i+1==st.dbidx[dbidx+1]) {
+            srmse_vec(dbidx)=std::sqrt(srmse/N);
+            srmse=0;
+            N=0;
+            dbidx++;
+        }
+    }
+    return srmse_vec;
+}
+
+t_type Analytics::calc_e_r_sq() const {
+    t_type e_r_sq_vec(st.dbidx.size()-1);
+    size_t dbidx=0;
+    t_type obs(st.size(dbidx));
+    t_type pred(st.size(dbidx));
+    size_t idx=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        obs(idx) = st(i).energy/st(i).natoms();
+        pred(idx) = stp(i).energy/st(i).natoms();
+        idx++;
+        if (i+1==st.dbidx[dbidx+1]) {
+            e_r_sq_vec(dbidx) = Statistics::r_sq(obs,pred);
+            dbidx++;
+            if (dbidx<st.dbidx.size()-2) {
+                obs.resize(st.size(dbidx));
+                pred.resize(st.size(dbidx));
+            }
+            idx=0;
+        }
+    }
+    return e_r_sq_vec;
+}
+
+t_type Analytics::calc_f_r_sq() const {
+    t_type f_r_sq_vec(st.dbidx.size()-1);
+    size_t dbidx=0;
+    t_type obs(3*st.calc_natoms(dbidx));
+    t_type pred(3*st.calc_natoms(dbidx));
+    size_t idx=0;
+    for (size_t i=0; i<st.size(); ++i) {
+        for (size_t a=0; a<st(i).natoms(); ++a) {
+            for (size_t k=0; k<3; ++k) {
+                obs(idx) = st(i).atoms[a].force(k);
+                pred(idx) = stp(i).atoms[a].force(k);
+                idx++;
+            }
+        }
+        if (i+1==st.dbidx[dbidx+1]) {
+            f_r_sq_vec(dbidx) = Statistics::r_sq(obs,pred);
+            dbidx++;
+            if (dbidx<st.dbidx.size()-2) {
+                obs.resize(3*st.calc_natoms(dbidx));
+                pred.resize(3*st.calc_natoms(dbidx));
+            }
+            idx=0;
+        }
+    }
+    return f_r_sq_vec;
+}
+
+t_type Analytics::calc_s_r_sq() const {
+    t_type s_r_sq_vec(st.dbidx.size()-1);
+    size_t dbidx=0;
+    t_type obs(6*st.size(dbidx));
+    t_type pred(6*st.size(dbidx));
+    size_t idx=0;
+    for (size_t i=0; i<st.size(); ++i) {
+
+        for (size_t x=0; x<3; ++x) {
+            for (size_t y=x; y<3; ++y) {
+                obs(idx) = st(i).stress(x,y);
+                pred(idx) = stp(i).stress(x,y);
+                idx++;
+            }
+        }
+
+        if (i+1==st.dbidx[dbidx+1]) {
+            s_r_sq_vec(dbidx) = Statistics::r_sq(obs,pred);
+            dbidx++;
+            if (dbidx<st.dbidx.size()-2) {
+                obs.resize(6*st.size(dbidx));
+                pred.resize(6*st.size(dbidx));
+            }
+            idx=0;
+        }
+    }
+    return s_r_sq_vec;
+}

analytics/analytics.h

+#ifndef ANALYTICS_H
+#define ANALYTICS_H
+
+#include "../structure_db.h"
+#include "../../CORE/core_types.h"
+
+
+/** Class for analysing and comparing datasets
+ *
+ */
+class Analytics {
+    private:
+        const StructureDB &st;    // orig data
+        const StructureDB &stp;   // predicted data
+
+    public:
+        Analytics(const StructureDB &st, const StructureDB &stp);
+
+        /** Return Energy/atom Mean Absolute Error for each DBFILE. */
+        t_type calc_e_mae() const;
+
+        /** Return Force Mean Absolute Error for each DBFILE. */
+        t_type calc_f_mae() const;
+
+        /** Return Stress Mean Absolute Error for each DBFILE.
+         *
+         *  Calculated using 6 independed components.
+         */
+        t_type calc_s_mae() const;
+
+        /** Return Energy/atom Root Mean Square Error for each DBFILE. */
+        t_type calc_e_rmse() const;
+
+        /** Return Force Root Mean Square Error for each DBFILE. */
+        t_type calc_f_rmse() const;
+
+        /** Return Stress Root Mean Square Error for each DBFILE.
+         *
+         *  Calculated using 6 independed components.
+         */
+        t_type calc_s_rmse() const;
+
+        /** Return Energy/atom coefficient of determination for each DBFILE. */
+        t_type calc_e_r_sq() const;
+
+        /** Return Force coefficient of determination for each DBFILE. */
+        t_type calc_f_r_sq() const;
+
+        /** Return Stress coefficient of determination for each DBFILE. */
+        t_type calc_s_r_sq() const;
+
+};
+#endif

analytics/statistics.cpp

+#include "statistics.h"
+#include <cmath>
+
+double Statistics::res_sum_sq(const vec &obs, const vec &pred) {
+    if (obs.size() != pred.size())
+        throw std::runtime_error("Containers differ in size.");
+
+    double rss=0;
+    for (size_t i=0; i<obs.size(); ++i) {
+        rss += std::pow(obs[i]-pred[i],2);
+    }
+    return rss;
+}
+
+double Statistics::tot_sum_sq(const vec &obs) {
+    double m = obs.mean();
+    double tss=0;
+    //for (auto v : obs) tss+=pow(v-m,2);
+    for (size_t i=0; i<obs.size(); ++i) tss+=pow(obs(i)-m,2);
+    return tss;
+
+}
+
+double Statistics::r_sq(const vec &obs, const vec &pred) {
+    return 1- res_sum_sq(obs,pred)/tot_sum_sq(obs);
+}
+
+double Statistics::variance(const vec &v) {
+    return tot_sum_sq(v)/(v.size()-1);
+}
+
+double Statistics::mean(const vec &v) {
+    return v.mean();
+}

analytics/statistics.h

+#ifndef STATISTICS_H
+#define STATISTICS_H
+
+#include "../../CORE/core_types.h"
+
+/** Some basis statistical tools */
+class Statistics {
+    using vec = aed_type2;
+    public:
+        /** Residual sum of squares. */
+        static double res_sum_sq(const vec &obs, const vec &pred);
+
+        /** Total sum of squares. */
+        static double tot_sum_sq(const vec &obs);
+
+        /** Coefficient of determination. */
+        static double r_sq(const vec &obs, const vec &pred);
+
+        /** Unbiased sample variance. */
+        static double variance(const vec &v);
+
+        /** Mean. */
+        static double mean(const vec &v);
+
+};
+#endif

atom.cpp

+#include "atom.h"
+#include "../CORE/periodic_table.h"
+#include <iomanip>
+
+Atom::Atom() {}
+
+Atom::Atom(const Element &element,
+        const double px, const double py, const double pz,
+        const double fx, const double fy, const double fz):
+    Element(element),
+    position(px,py,pz),
+    force(fx,fy,fz)
+{}
+
+std::ostream& operator<<(std::ostream& os, const Atom& atom)
+{
+    os     << (Element&)atom;
+    os     << "Coordinates:  " << std::left << atom.position << std::endl;
+    os     << "Force:        " << std::left << atom.force << std::endl;
+    return os;
+}
+bool Atom::operator==(const Atom &a) const {
+    double EPSILON = std::numeric_limits<double>::epsilon();
+    return 
+        position.isApprox(a.position, EPSILON)
+        && force.isApprox(a.force, EPSILON)
+        && Element::operator==(a)
+        ;
+}

atom.h

+#ifndef ATOM_H
+#define ATOM_H
+
+#include <string>
+#include <iostream>
+#include "../CORE/element.h"
+#include "../CORE/core_types.h"
+
+/**
+ * Container to represent atom properties
+ *
+ * Usage example:
+ *
+ *     # create an empty atom object - all attributes are left uninitialised
+ *     Atom atom;
+ *
+ *     # set atom symbol
+ *     atom.symbol="Ti"
+ *
+ *     # set atom name, atomic number... see Element class for attributes
+ *
+ *     # set atom position to (1.1, 2.2, 3.3)
+ *     atom.position(0) = 1.1;
+ *     atom.position(1) = 2.2;
+ *     atom.position(2) = 3.3;
+ *
+ *     # set force
+ *     atom.force(0)= 0.1;
+ *     atom.force(1)= -0.2;
+ *     atom.force(2)= 0.3;
+ *
+ * Usage example:
+ *
+ *     # Use constructor to fully initialise this object
+ *     # with the position and force as in the example above
+ *     Element element = PeriodicTable().find_by_symbol("Ti");
+ *     Atom atom(element, 1.1, 2.2, 3.3, 0.1, -0.2, 0.3);
+ *
+ * Usage example:
+ *
+ *     # Print atom object using streams:
+ *     std::cout << atom;
+ *
+ *     # Print position only
+ *     std::cout << atom.position
+ *
+ * @see Structure
+ */
+struct Atom: public Element {
+    /**
+     * Create an empty atom object. All class attributes are left uninitialised.
+     */
+    Atom();
+
+    /** This constructor fully initialise this object
+     *
+     * @param[in] element   Chemical Element
+     * @param[in] px,py,pz  Atomic coordinates
+     * @param[in] fx,fy,fz  Force acting on the atom
+     *
+     */
+    Atom(const Element &element,
+            const double px, const double py, const double pz,
+            const double fx, const double fy, const double fz);
+
+    /** Hold position of the atom. */
+    Vec3d position;
+
+    /** Hold force on the atom. */
+    Vec3d force;
+
+    /** Print object summary to the stream */
+    friend std::ostream& operator<<(std::ostream& os, const Atom& atom);
+
+    /** Return true if both atoms are the same.
+     *
+     * This operator compares chemical type, position and force.
+     */
+    bool operator==(const Atom &) const;
+
+};
+#endif

descriptors_calc.h

+#ifndef DESCRIPTORS_CALC_H
+#define DESCRIPTORS_CALC_H
+
+#include "descriptors_calc_base.h"
+#include "../CORE/config/config.h"
+#include "structure.h"
+#include "structure_db.h"
+#include "st_descriptors.h"
+#include "st_descriptors_db.h"
+#include "../MODELS/descriptors/d2/d2_base.h"
+#include "../MODELS/descriptors/d3/d3_base.h"
+#include "../MODELS/descriptors/dm/dm_base.h"
+#include "../MODELS/cutoffs/cutoffs.h"
+
+/** \brief Descriptors calculator
+ *
+ * This object does not store any data.
+ *
+ * Usage example:
+ * \code{.cpp}
+ * # Build calculator using Config object.
+ * DescriptorsCalc<D2_LJ, D3_Dummy, DM_Dummy> dc(config);
+ *
+ * # Calculate all descriptors in a Structure st
+ * StDescriptors std = dc.calc(st);
+ *
+ * # Calculate all descriptors in a StructureDB st_db
+ * StDescriptorsDB std_db = dc.calc(st_db);
+ * \endcode
+ *
+ *  This object sets the following \ref INTERNAL_KEYS
+ *  to the provided Config object:
+ *
+ *  \ref SIZE2B,\ref SIZE3B,\ref SIZEMB,\ref DSIZE
+ *
+ *  Required keys:
+ *  \ref INIT2B,\ref INIT3B,\ref INITMB,
+ *  \ref RCUT2B,\ref RCUT3B,\ref RCUTMB (for initialised ones)
+ *
+ *
+ *  @tparam D2 D2_Base child, two-body type calculator.
+ *  @tparam D3 D3_Base child, three-body type calculator.
+ *  @tparam DM DM_Base child, many-body type calculator.
+ *  @tparam C2 Cut_Base child, two-body cutoff function.
+ *  @tparam C3 Cut_Base child, three-body cutoff function.
+ *  @tparam CM Cut_Base child, many-body cutoff function.
+ */
+template <typename D2=D2_Base&, typename D3=D3_Base&, typename DM=DM_Base&, typename C2=Cut_Base&, typename C3=Cut_Base&, typename CM=Cut_Base&>
+class DescriptorsCalc: public DC_Base {
+
+    public:
+        Config &config;
+        /** Constructor to fully initialise this object.
+         *
+         *  REQUIRED CONFIG KEYS:
+         *  \ref FORCE, \ref STRESS,
+         *  at least one of: \ref RCUT2B, \ref RCUT3B,\ref RCUTMB
+         *
+         */
+        DescriptorsCalc(Config &c);
+
+        /** This constructor is equivalent to the main constructor above.
+         *
+         * The difference is technical in nature.
+         * Here we use OOP with generic programming to call virtual functions
+         * to calculate descriptors - this incurs small efficiency penalty
+         * but makes selecting of descriptors somewhat easier at runtime.
+         *
+         * In contrast, the main constructor above is purely generic.
+         *
+         * The rule of thumb is: Use the main contructor unless it's a pain.
+         *
+         * Usage example:
+         * \code{.cpp}
+         * # Build calculator using Config object.
+         * D2_Base *d2 = new D2_LJ(config);
+         * D2_Base *d3 = new D3_Dummy(config);
+         * D2_Base *dm = new DM_Dummy(config);
+         * DescriptorsCalc<> dc(config, d2, d3, dm);
+         *
+         * # Calculate all descriptors in a Structure st
+         * StDescriptors std = dc.calc(st);
+         *
+         * # Calculate all descriptors in a StructureDB st_db
+         * StDescriptorsDB std_db = dc.calc(st_db);
+         * \endcode
+         *
+         * \warning {
+         *  The unfortunate consequence is that
+         *  you might want to try this
+         *  \code{.cpp}
+         *  DescriptorsCalc<> dc(config);
+         *  \endcode
+         *  Basically you are trying to init abstract classes with config object.
+         *  It won't compile!
+         * }
+         *
+         */
+        template <typename T1, typename T2, typename T3,typename T4, typename T5, typename T6>
+        DescriptorsCalc(Config &c, T1 &d2, T2 &d3, T3 &dm, T4 &c2, T5 &c3, T6 &cm);
+
+        /** Calculate all descriptors in a Structure st */
+        StDescriptors calc(const Structure &st);
+
+        /** Calculate all descriptors in a StructureDB st_db */
+        StDescriptorsDB calc(const StructureDB &st_db);
+
+
+        /** Calculate density vector for the structure */
+        void calc_rho(const Structure &st, StDescriptors &std);
+
+        C2 c2;
+        C3 c3;
+        CM cm;
+        D2 d2;
+        D3 d3;
+        DM dm;
+    private:
+        template <typename T1, typename T2, typename T3>
+        DescriptorsCalc(Config &c,T1 &t1,T2 &t2,T3 &t3);
+
+        void calc(const Structure &st, StDescriptors &std);
+        void calc_dimer(const Structure &st, StDescriptors &std);
+        void common_constructor();
+
+};
+
+#include "descriptors_calc.hpp"
+#endif

descriptors_calc.hpp

+#ifndef DESCRIPTORS_CALC_HPP
+#define DESCRIPTORS_CALC_HPP
+
+// This file is includeed back to descriptors_calc.h
+// So There is no need to include headers
+// Useful for debugging though...
+#include "descriptors_calc.h"
+
+#include <cstdio>
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+template <typename T1, typename T2, typename T3>
+DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c, T1 &t1, T2 &t2, T3 &t3):
+    config(c),
+    d2(t1),
+    d3(t2),
+    dm(t3)
+{
+    if (!config.exist("DSIZE"))
+        common_constructor();
+    else {
+        size_t bias=0;
+        if (config.get<bool>("BIAS"))
+            bias++;
+        if (config.get<bool>("INIT2B")) {
+            d2.fidx=bias;
+        }
+        if (config.get<bool>("INITMB")) {
+            dm.fidx=bias+d2.size();
+        }
+    }
+}
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c):
+    DescriptorsCalc(c,c,c,c)
+{
+    if (c.get<bool>("INIT2B")) {
+        c2 = C2(c.get<double>("RCUT2B"));
+        if (!config.exist("RCTYPE2B"))
+            config.add("RCTYPE2B",c2.label());
+    }
+    else {
+        c2 = C2(0);
+    }
+    if (c.get<bool>("INIT3B")) {
+        c3 = C3(c.get<double>("RCUT3B"));
+        if (!config.exist("RCTYPE3B"))
+            config.add("RCTYPE3B",c3.label());
+    }
+    else {
+        c3 = C3(0);
+    }
+    if (c.get<bool>("INITMB")) {
+        cm = CM(c.get<double>("RCUTMB"));
+        if (!config.exist("RCTYPEMB"))
+            config.add("RCTYPEMB",cm.label());
+    }
+    else {
+        cm = CM(0);
+    }
+}
+
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+        template <typename T1, typename T2, typename T3,typename T4, typename T5, typename T6>
+DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c, T1 &d2, T2 &d3, T3 &dm, T4 &c2, T5 &c3, T6 &cm):
+    config(c),
+    c2(c2),
+    c3(c3),
+    cm(cm),
+    d2(d2),
+    d3(d3),
+    dm(dm)
+{
+    if (c.get<bool>("INIT2B")) {
+        if (!config.exist("RCTYPE2B"))
+            config.add("RCTYPE2B",c2.label());
+    }
+    if (c.get<bool>("INIT3B")) {
+        if (!config.exist("RCTYPE3B"))
+            config.add("RCTYPE3B",c3.label());
+    }
+    if (c.get<bool>("INITMB")) {
+        if (!config.exist("RCTYPEMB"))
+            config.add("RCTYPEMB",cm.label());
+    }
+
+    if (!config.exist("DSIZE"))
+        common_constructor();
+}
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::common_constructor() {
+    size_t dsize=0;
+
+    // Calculate total size of the descriptor.
+    // Add relevant keys to the config.
+    // Set first indices for all types
+    size_t bias=0;
+    if (config.get<bool>("BIAS"))
+        bias++;
+    if (config.get<bool>("INIT2B")) {
+        config.add("SIZE2B",d2.size());
+        if (!config.exist("TYPE2B"))
+            config.add("TYPE2B",d2.label());
+        dsize+=d2.size();
+        d2.fidx=bias;
+    }
+    else {
+        config.add("SIZE2B",0);
+    }
+
+    if (config.get<bool>("INIT3B")) {
+        config.add("SIZE3B",d3.size());
+        if (!config.exist("TYPE3B"))
+            config.add("TYPE3B",d3.label());
+        dsize+=d3.size();
+    }
+    else {
+        config.add("SIZE3B",0);
+    }
+
+    if (config.get<bool>("INITMB")) {
+        config.add("SIZEMB",dm.size());
+        if (!config.exist("TYPEMB"))
+            config.add("TYPEMB",dm.label());
+        dsize+=dm.size();
+        dm.fidx=bias+d2.size();
+    }
+    else {
+        config.add("SIZEMB",0);
+    }
+
+    if (!dsize)
+        throw std::runtime_error("The descriptor size is 0, check your config.");
+
+    if (config.get<bool>("BIAS"))
+        dsize++;
+
+    config.add("DSIZE",dsize);
+}
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_rho(const Structure &st, StDescriptors &st_d) {
+    double rcut_mb_sq = pow(config.get<double>("RCUTMB"),2);
+    rhos_type &rhos = st_d.rhos;
+    size_t s = dm.rhoi_size()+dm.rhoip_size();
+    rhos.resize(s,st.natoms());
+    rhos.set_zero();
+
+  Vec3d delij;
+    for (size_t i=0; i<st.natoms(); ++i) {
+        const Atom &a1 = st(i);
+
+        for (size_t jj=0; jj<st.nn_size(i); ++jj) {
+            const Vec3d &a2pos = st.nn_pos(i,jj);
+            //delij = a1.position - a2pos;
+            delij[0] = a1.position[0] - a2pos[0];
+            delij[1] = a1.position[1] - a2pos[1];
+            delij[2] = a1.position[2] - a2pos[2];
+
+            //double rij_sq = delij * delij;
+            double rij_sq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
+            if (rij_sq > rcut_mb_sq) continue;
+            int Zj = st.near_neigh_atoms[i][jj].Z;
+            double rij = sqrt(rij_sq);
+            double fc_ij = Zj*cm.calc(rij);
+            dm.calc_rho(rij,rij_sq,fc_ij,delij,st_d.get_rho(i));
+        }
+    }
+}
+
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, StDescriptors &st_d) {
+
+    bool init2b = config.get<bool>("INIT2B");
+    bool initmb = config.get<bool>("INITMB");
+
+    size_t size2b=0;
+    size_t sizemb=0;
+    if (init2b) size2b=config.get<size_t>("SIZE2B");
+    if (initmb) sizemb=config.get<size_t>("SIZEMB");
+
+    // deal with the case initmb is set to true but dummy is used
+    init2b = init2b && size2b;
+    initmb = initmb && sizemb;
+
+    size_t bias=0;
+    if (config.get<bool>("BIAS"))
+        bias++;
+
+    if (initmb)
+        calc_rho(st,st_d);
+
+    double rcut_max_sq = pow(config.get<double>("RCUTMAX"),2);
+    double rcut_2b_sq = 0.0;
+    double rcut_mb_sq = 0.0;
+
+    if (init2b) rcut_2b_sq = pow(config.get<double>("RCUT2B"),2);
+    if (initmb) rcut_mb_sq = pow(config.get<double>("RCUTMB"),2);
+
+    // zero all aeds and set bias
+    for (size_t i=0; i<st.natoms(); ++i) {
+        aed_type2 &aed = st_d.get_aed(i);
+        aed.set_zero();
+        aed(0)=static_cast<double>(bias);   // set bias
+    }
+
+    // calculate many-body energy
+    // do it before main loop so rho prime
+    // can be calculated
+    if (initmb) {
+        for (size_t i=0; i<st.natoms(); ++i) {
+            aed_type2 &aed = st_d.get_aed(i);
+            rho_type& rhoi = st_d.get_rho(i);
+            dm.calc_aed(rhoi,aed);
+        }
+    }
+
+  bool use_force = config.get<bool>("FORCE");
+  bool use_stress = config.get<bool>("STRESS");
+
+  Vec3d delij;
+    for (size_t i=0; i<st.natoms(); ++i) {
+        const Atom &a1 = st(i);
+        aed_type2 &aed = st_d.get_aed(i);
+
+        for (size_t jj=0; jj<st.nn_size(i); ++jj) {
+            const Vec3d &a2pos = st.nn_pos(i,jj);
+
+            delij[0] = a1.position[0] - a2pos[0];
+            delij[1] = a1.position[1] - a2pos[1];
+            delij[2] = a1.position[2] - a2pos[2];
+
+            double rij_sq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
+
+            if (rij_sq > rcut_max_sq) continue;
+            int Zj = st.near_neigh_atoms[i][jj].Z;
+            double rij = sqrt(rij_sq);
+            double rij_inv = 1.0/rij;
+
+            // CALCULATE TWO-BODY TERM
+            if (use_force || use_stress) {
+                fd_type &fd_ij = st_d.fd[i][jj];
+                if (rij_sq <= rcut_2b_sq && init2b) {
+                    double fc_ij = Zj*c2.calc(rij);
+                    double fcp_ij = Zj*c2.calc_prime(rij);
+                    d2.calc_all(rij,rij_sq,fc_ij,fcp_ij,aed,fd_ij);
+                    // Two-body descriptor calculates x-direction only - fd_ij(n,0)
+                    // so we have to copy x-dir to y- and z-dir
+                    // and scale them by the unit directional vector delij/rij.
+                    for (size_t n=bias; n<size2b+bias; ++n) {
+                        fd_ij(n,0) *= 0.5*rij_inv;
+                        fd_ij(n,1) = fd_ij(n,0)*delij[1];
+                        fd_ij(n,2) = fd_ij(n,0)*delij[2];
+                        fd_ij(n,0) *= delij[0];
+                    }
+                }
+                // CALCULATE MANY-BODY TERM
+                if (rij_sq <= rcut_mb_sq && initmb) {
+                    double fc_ij = Zj*cm.calc(rij);
+                    double fcp_ij = Zj*cm.calc_prime(rij);
+                    rho_type& rhoi = st_d.get_rho(i);
+                    int mode = dm.calc_dXijdri(rij,rij_sq,delij,
+                            fc_ij,fcp_ij,rhoi,fd_ij);
+                    if (mode==0) {
+                        // some dm compute x-dir only, similarly to d2 above
+                        for (size_t n=size2b+bias; n<size2b+sizemb+bias; ++n) {
+                            fd_ij(n,0) *= rij_inv;
+                            fd_ij(n,1) = fd_ij(n,0);
+                            fd_ij(n,2) = fd_ij(n,0);
+                            fd_ij(n,0) *= delij[0];
+                            fd_ij(n,1) *= delij[1];
+                            fd_ij(n,2) *= delij[2];
+                        }
+                    }
+                }
+            }
+            else {
+                if (rij_sq <= rcut_2b_sq && init2b) {
+                    double fc_ij = Zj*c2.calc(rij);
+                    d2.calc_aed(rij,rij_sq,fc_ij,aed);
+                }
+            }
+
+        }
+    }
+    if (init2b) {
+        for (size_t n=0; n<st.natoms(); ++n) {
+            for(size_t s=bias; s<bias+d2.size(); ++s) {
+                st_d.get_aed(n)(s) *= 0.5;
+            }
+        }
+    }
+}
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescriptors &st_d) {
+
+    double r_b = config.get<double>("DIMER",1);
+    bool bond_bool = config.get<bool>("DIMER",2);
+
+    // Here we assume that training data consists of 4 atoms
+    // The cutoff distance specified in the config file works in the usual way
+    // i.e. the maximum distance between two atoms 
+    // I,J are molecules, i1,i2 atom in I, similarly for J
+    //    I- - - - - - - - - - - -J
+    // i1---i2- - - - - - - - -j1---j2
+    //    |- - - -r_com- - - - - -|   
+    //  |- - - - -r_max- - - - - - -|
+
+    bool init2b = config.get<bool>("INIT2B");
+    bool initmb = config.get<bool>("INITMB");
+    size_t size2b=0;
+    size_t sizemb=0;
+    if (init2b) size2b=config.get<size_t>("SIZE2B");
+    if (initmb) sizemb=config.get<size_t>("SIZEMB");
+
+    // deal with the case initmb is set to true and dummy is used
+    init2b = init2b && size2b;
+    initmb = initmb && sizemb;
+
+    //double rcut_max_sq = pow(config.get<double>("RCUTMAX"),2);
+    double rcut_2b_sq = 0.0;
+    double rcut_mb_sq = 0.0;
+
+    if (init2b) rcut_2b_sq = pow(config.get<double>("RCUT2B"),2);
+    if (initmb) rcut_mb_sq = pow(config.get<double>("RCUTMB"),2);
+
+    // Max distance between CoM of two interacting molecules
+    double rcut_com_sq = pow(config.get<double>("RCUTMAX")-r_b,2);
+
+    // Not that this differ to how lammps implements this
+    // TODO make it consistent between those two
+    Mat6R3C delM;  //i1-j1,i2-j1,i1-j2,i2-j2,i1-i2,j1-j2
+    double r_sq[6];
+    double r[6];
+
+    const std::vector<Atom> &atoms = st.atoms;
+
+    size_t bias=0;
+    if (config.get<bool>("BIAS"))
+        bias++;
+
+    // TODO weighting factors
+    // For now assume all are the same type
+    //int Zj = st.near_neigh_atoms[0][0].Z;
+    int Zj = 1;
+
+    // map of atom label and distances
+    // idx[0] - distance label between 0 and 2 atom
+    std::vector<std::pair<int,int>> idx(6);
+    idx[0] = std::make_pair(0,2);
+    idx[1] = std::make_pair(1,2);
+    idx[2] = std::make_pair(0,3);
+    idx[3] = std::make_pair(1,3);
+    idx[4] = std::make_pair(0,1);
+    idx[5] = std::make_pair(2,3);
+
+    // zero all aeds+rho and set bias
+    for (size_t i=0; i<4; ++i) {
+        aed_type2 &aed = st_d.get_aed(i);
+        aed.set_zero();
+        aed(0)=static_cast<double>(bias);   // set bias
+    }
+    if (initmb) {
+        size_t s = dm.rhoi_size()+dm.rhoip_size();
+        st_d.rhos.resize(s,4);
+        st_d.rhos.set_zero();
+    }
+
+    Vec3d xicom = 0.5*(atoms[0].position + atoms[1].position);
+    Vec3d xjcom = 0.5*(atoms[2].position + atoms[3].position);
+
+    Vec3d del_com = xicom-xjcom;;
+    double r_com_sq = del_com * del_com;
+    if (r_com_sq > rcut_com_sq) {
+        return;
+    }
+
+    // compute all 6 distances
+    for (size_t n=0; n<6; ++n) {
+        delM.row(n) = atoms[idx[n].first].position - atoms[idx[n].second].position;
+        r_sq[n] =  delM.row(n) * delM.row(n);
+        r[n] = sqrt(r_sq[n]);
+    }
+
+    // compute densities and 2b aed for every atom
+    // if bond is not included use 0,1,2,3 distances,
+    // otherwise use all
+    size_t N = bond_bool ? 6 : 4;
+    for (size_t n=0; n<N; ++n) {
+        if (r_sq[n] <= rcut_mb_sq && initmb) {
+            double fcmb = Zj*cm.calc(r[n]);
+            dm.calc_rho(r[n],r_sq[n],fcmb,delM.row(n),st_d.get_rho(idx[n].first));
+            dm.calc_rho(r[n],r_sq[n],fcmb,-delM.row(n),st_d.get_rho(idx[n].second));
+        }
+        // Do not compute 2b term between bonded atoms
+        if (r_sq[n] <= rcut_2b_sq && init2b) {
+            double fc2b = Zj*c2.calc(r[n]);
+            d2.calc_aed(r[n],r_sq[n],fc2b,st_d.get_aed(idx[n].first));
+            d2.calc_aed(r[n],r_sq[n],fc2b,st_d.get_aed(idx[n].second));
+        }
+    }
+
+    // calculate many-body aed
+    if (initmb) {
+        for (size_t i=0; i<4; ++i) {
+            aed_type2 &aed = st_d.get_aed(i);
+            rho_type& rhoi = st_d.get_rho(i);
+            dm.calc_aed(rhoi,aed);
+        }
+    }
+
+    //if (init2b) {
+    //    for (size_t n=0; n<st.natoms(); ++n) {
+    //        for(size_t s=bias; s<bias+d2.size(); ++s) {
+    //            aeds[n](s) *= 0.5;
+    //        }
+    //    }
+    //}
+}
+
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+StDescriptors DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st) {
+    StDescriptors st_d(st, config);
+    if (config.get<bool>("DIMER",0))
+        calc_dimer(st,st_d);
+    else
+        calc(st,st_d);
+    return st_d;
+}
+
+template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
+StDescriptorsDB DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const StructureDB &stdb) {
+    StDescriptorsDB st_desc_db(stdb, config);
+
+#pragma omp parallel for
+    for(size_t i=0; i<stdb.size(); ++i)
+    if (config.get<bool>("DIMER",0))
+        calc_dimer(stdb(i),st_desc_db(i));
+    else
+        calc(stdb(i),st_desc_db(i));
+    return st_desc_db;
+}
+#endif

descriptors_calc_base.h

+#ifndef DC_Base_H
+#define DC_Base_H
+
+#include "structure.h"
+#include "structure_db.h"
+#include "st_descriptors.h"
+#include "st_descriptors_db.h"
+
+class DC_Base {
+    public:
+        virtual ~DC_Base() {};
+        virtual StDescriptors calc(const Structure &) { return StDescriptors(); };
+        virtual StDescriptorsDB calc(const StructureDB &) {return StDescriptorsDB();};
+};
+#endif

design_matrix/design_matrix.h

+#ifndef DESIGN_MATRIX_H
+#define DESIGN_MATRIX_H
+
+#include <stdexcept>
+#include "../../CORE/config/config.h"
+#include "../st_descriptors_db.h"
+#include "../structure_db.h"
+#include "../descriptors_calc.h"
+#include "../normaliser.h"
+
+
+/** Design Matrix \f$ \Phi \f$ and the corresponding target vector \f$ \mathrm{T} \f$
+ *
+ * An aggregation matrix composed of \phi_i matrices
+ *
+ *     Phi = | \phi_1 |
+ *           | \phi_2 |
+ *           | ...    |
+ *           | \phi_N |
+ *
+ * where N corresponds to the number of structures in a provided StructureDB.
+ *
+ * The ordering in \f$ \phi_i \f$
+ *
+ *     | energy of the 1st structure
+ *     | force on the 1st atom in the x-dir
+ *     | force on the 1st atom in the y-dir
+ *     | force on the 1st atom in the z-dir
+ *     | force on the 2st atom in the x-dir
+ *     | ...
+ *     | force on the n-th atom in the z-dir
+ *     | stress xx
+ *     | stress xy
+ *     | stress xz
+ *     | stress yy
+ *     | stress yz
+ *     | stress zz
+ *
+ * We need to know: \ref STRESS, \ref FORCE
+ *
+ * Scalling factors \ref FWEIGHT, \ref  SWEIGHT
+ *
+ * Default scalling:
+ *  - Forces are scaled by 1/N_i/3 where N_i is a number of atoms in the i-th structure.
+ *  - Stresses are scaled by 1/6.
+ *  - Energies are scaled by 1/N_i
+ *
+ * @tparam F Function_Base child -> BF_Base or Kern_Base
+ *
+ */
+template <typename F>
+class DesignMatrix {
+
+    public:
+
+        F f;
+        phi_type Phi;
+        t_type T;
+        t_type Tlabels; // 0-Energy, 1-Force, 2-Stress
+        bool scale=true;    // Control escale,fscale,sscale
+
+        double e_std_dev=1;
+        double f_std_dev=1;
+        double s_std_dev[6] = {1,1,1,1,1,1};
+
+        /** \brief This constructor fully initialise this object
+         *
+         * This class is used to build a Design Matrix.
+         *
+         * Usage example:
+         * \code{.cpp}
+         * Config config("Config");
+         * BF_Linear bf(config);
+         * DesignMatrix<LinearKernel> desmat(bf, config);
+         *
+         * \endcode
+         */
+        DesignMatrix(F &f, Config &c):
+            f(f),
+            config(c),
+            force(config.get<bool>("FORCE")),
+            stress(config.get<bool>("STRESS")),
+            verbose(config.get<int>("VERBOSE")),
+            eweightglob(config.get<double>("EWEIGHT")),
+            fweightglob(config.get<double>("FWEIGHT")),
+            sweightglob(config.get<double>("SWEIGHT"))
+    {
+        if (config.exist("ESTDEV"))
+            e_std_dev = config.get<double>("ESTDEV");
+        if (config.exist("FSTDEV"))
+            f_std_dev = config.get<double>("FSTDEV");
+        if (config.exist("SSTDEV"))
+            config.get<double[6]>("SSTDEV", s_std_dev);
+    }
+
+
+        /** \brief Build design matrix from already calculated StDescriptorsDB
+         *
+         * Here we simply build matrix from already calculated descriptors.
+         * The vector of targets  **T** is build from StructureDB.
+         * D2, D3 and DM calculators are not used.
+         */
+        void build(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
+            resize(stdb);
+            compute_stdevs(stdb);
+            fill_T(stdb);
+            std::vector<size_t> rows(stdb.size());
+            size_t row=0;
+            for (size_t s=0; s<stdb.size(); ++s) {
+                rows[s] = row;
+                row++;
+                if (force) {
+                    for (size_t a=0; a<stdb(s).natoms(); ++a) {
+                        row+=3;
+                    }
+                }
+                if (stress) {
+                    row+=6;
+                }
+            }
+            // loop over all structures and build
+#pragma omp parallel for
+            for (size_t s=0; s<stdb.size(); ++s) {
+                build(rows[s], stdb(s),st_desc_db(s));
+            }
+        };
+
+        /** \brief Calculate descriptors and build design matrix. */
+        template <typename DC>
+        void build(const StructureDB &stdb, Normaliser &norm,
+                DC &dc) {
+            //DescriptorsCalc<D2,D3,DM,C2,C3,CM> dc(config);
+            resize(stdb);   // call after dc set DSIZE key
+            compute_stdevs(stdb);
+            fill_T(stdb);
+            // for opm we need to find first rows for each structure
+
+            std::vector<size_t> rows(stdb.size());
+            size_t row=0;
+            for (size_t s=0; s<stdb.size(); ++s) {
+                rows[s] = row;
+                row++;
+                if (force) {
+                    for (size_t a=0; a<stdb(s).natoms(); ++a) {
+                        row+=3;
+                    }
+                }
+                if (stress) {
+                    row+=6;
+                }
+            }
+#pragma omp parallel for
+            for (size_t s=0; s<stdb.size(); ++s) {
+                StDescriptors st_d = dc.calc(stdb(s));
+
+                if(config.get<bool>("NORM"))
+                    norm.normalise(st_d);
+
+                build(rows[s], stdb(s), st_d);
+            }
+
+        }
+
+    private:
+        Config &config;
+        size_t rows=0;
+        size_t cols=0;
+        bool force;
+        bool stress;
+        int verbose;
+        double eweightglob;
+        double fweightglob;
+        double sweightglob;
+
+        // resize Phi and T and set all elements to zero
+        void resize(const StructureDB &stdb) {
+            rows=0;
+            cols = f.get_phi_cols(config);
+
+            for (size_t s=0; s<stdb.size(); ++s) {
+                rows++;
+                if (stress) rows+=6;
+                if (force) rows+=stdb(s).natoms()*3;
+            }
+
+            T.resize(rows);
+            Tlabels.resize(rows);
+            Phi.resize(rows,cols);
+
+            if (verbose) {
+                std::cout << "Phi rows: "<< Phi.rows() << std::endl;
+                std::cout << "Phi cols: "<< Phi.cols() << std::endl;
+            }
+
+            T.set_zero();
+            Phi.set_zero();
+
+        }
+
+        void compute_stdevs(const StructureDB &stdb) {
+            t_type evec(stdb.size());
+            Matrix svec(stdb.size(),6);
+            size_t num_forces=0;
+            for (size_t s=0; s<stdb.size(); ++s) {
+                evec(s) = stdb(s).energy/stdb(s).natoms();
+                if (stress) {
+                    size_t j=0;
+                    for (size_t x=0; x<3; ++x) {
+                        for (size_t y=x; y<3; ++y) {
+                            svec(s,j++)=stdb(s).stress(x,y);
+                        }
+                    }
+                }
+                num_forces +=3*stdb(s).natoms();
+            }
+
+            //e_std_dev = std::sqrt((evec - evec.mean()).square().sum()/(evec.size()-1));
+            e_std_dev = evec.std_dev(evec.mean(), evec.size()-1);
+            if (verbose) std::cout << "Energy standard deviation (eV/atom): " << e_std_dev << std::endl;
+
+            if (stress) {
+                for (size_t j=0; j<6; ++j) {
+                    s_std_dev[j] = svec.col(j).std_dev(svec.col(j).mean(), svec.col(j).size()-1);
+                    if (verbose) std::cout << "Stress standard deviation (eV): " << s_std_dev[j] << std::endl;
+                }
+            }
+
+
+            if (force) {
+                size_t j=0;
+                t_type fvec(num_forces);
+                for (size_t s=0; s<stdb.size(); ++s) {
+                    for (const Atom &a : stdb(s).atoms) {
+                        fvec(j++) = a.force(0);
+                        fvec(j++) = a.force(1);
+                        fvec(j++) = a.force(2);
+                    }
+                }
+                //fvec /= e_std_dev;
+                //svec /= e_std_dev;
+                // e_std_dev has units of energy
+                // f_std_dev has units of inverse distance
+                f_std_dev = fvec.std_dev(fvec.mean(),fvec.size()-1);
+                if (verbose) std::cout << "Force standard deviation (1/A): " << f_std_dev << std::endl;
+            }
+
+            config.add("ESTDEV",e_std_dev);
+            config.add("FSTDEV",f_std_dev);
+            for (size_t j=0; j<6; ++j)
+                config.add("SSTDEV",s_std_dev[j]);
+        }
+        void fill_T(const StructureDB &stdb) {
+            size_t j=0;
+            for (size_t s=0; s<stdb.size(); ++s) {
+
+
+                double escale = 1;
+
+                if (scale) escale = stdb(s).eweight*eweightglob/stdb(s).natoms();///e_std_dev;
+                Tlabels(j) = 0;
+                T(j++) = stdb(s).energy*escale;
+
+                if (force) {
+                    double fscale = 1;
+                    if (scale) fscale = stdb(s).fweight*fweightglob/stdb(s).natoms()/3.0;///e_std_dev/f_std_dev;
+                    for (const Atom &a : stdb(s).atoms) {
+                        Tlabels(j) = 1;
+                        T(j++) = a.force(0)*fscale;
+                        Tlabels(j) = 1;
+                        T(j++) = a.force(1)*fscale;
+                        Tlabels(j) = 1;
+                        T(j++) = a.force(2)*fscale;
+                    }
+                }
+                if (stress) {
+                    double sscale_arr[6] {1,1,1,1,1,1};
+                    if (scale)
+                        for(size_t xy=0;xy<6;++xy)
+                            sscale_arr[xy] = stdb(s).sweight*sweightglob/6.0;///e_std_dev/s_std_dev[xy];
+                    size_t xy=0;
+                    for (size_t x=0; x<3; ++x) {
+                        for (size_t y=x; y<3; ++y) {
+                            Tlabels(j) = 2;
+                            T(j++)=stdb(s).stress(x,y)*sscale_arr[xy++];
+                        }
+                    }
+                }
+            }
+        }
+        void build(size_t &row, const Structure &st, const StDescriptors &st_d) {
+
+            double escale = 1;
+
+            if (scale) escale = st.eweight*eweightglob/st.natoms();
+            f.calc_phi_energy_row(Phi,row,escale,st,st_d);
+            if (force) {
+                double fscale = 1;
+                if (scale) fscale = st.fweight*fweightglob/st.natoms()/3.0;///f_std_dev;
+                f.calc_phi_force_rows(Phi,row,fscale,st,st_d);
+            }
+            if (stress) {
+                double sscale_arr[6] {1,1,1,1,1,1};
+                if (scale)
+                    for(size_t xy=0;xy<6;++xy)
+                        sscale_arr[xy] = st.sweight*sweightglob/6.0;///s_std_dev[xy];
+                f.calc_phi_stress_rows(Phi,row,sscale_arr,st,st_d);
+            }
+        }
+};
+#endif

design_matrix/functions/basis_functions/dm_bf_all.h

+#include "dm_bf_linear.h"
+#include "dm_bf_polynomial2.h"

design_matrix/functions/basis_functions/dm_bf_base.cpp

+#include "dm_bf_base.h"
+DM_BF_Base::~DM_BF_Base() {}

design_matrix/functions/basis_functions/dm_bf_base.h

+#ifndef DM_BASIS_FUNCTIONS_H
+#define DM_BASIS_FUNCTIONS_H
+
+#include <iostream>
+#include "../dm_function_base.h"
+#include "../../../../CORE/core_types.h"
+#include "../../../../MODELS/functions/basis_functions/bf_base.h"
+#include "../../../structure.h"
+#include "../../../st_descriptors.h"
+
+struct DM_BF_Base: public DM_Function_Base, public virtual BF_Base {
+    virtual ~DM_BF_Base();
+    virtual size_t get_phi_cols(const Config &config)=0;
+    virtual void calc_phi_energy_row(phi_type &Phi, size_t &row,
+            const double fac, const Structure &st, const StDescriptors &st_d)=0;
+    virtual void calc_phi_force_rows(phi_type &Phi, size_t &row,
+            const double fac, const Structure &st, const StDescriptors &st_d)=0;
+    virtual void calc_phi_stress_rows(phi_type &Phi, size_t &row,
+            const double fac[6], const Structure &st, const StDescriptors &st_d)=0;
+};
+#endif

design_matrix/functions/basis_functions/dm_bf_linear.cpp

+#include "dm_bf_linear.h"
+
+Registry<DM_Function_Base>::Register<DM_BF_Linear> DM_BF_Linear_1( "BF_Linear" );
+Registry<DM_Function_Base,Config&>::Register<DM_BF_Linear> DM_BF_Linear_2( "BF_Linear" );
+
+DM_BF_Linear::DM_BF_Linear() {}
+DM_BF_Linear::DM_BF_Linear(const Config &c): BF_Linear(c)
+{}
+size_t DM_BF_Linear::get_phi_cols(const Config &config)
+{
+  size_t cols = config.get<size_t>("DSIZE");
+  return cols;
+}
+void DM_BF_Linear::calc_phi_energy_row(phi_type &Phi, size_t &row,
+                                       const double fac, const Structure &, const StDescriptors &st_d)
+{
+  for (size_t i=0; i<st_d.naed(); ++i) {
+    const aed_type2 &aed = st_d.get_aed(i);
+    for (size_t j=0; j<aed.size(); ++j) {
+      Phi(row,j)+=aed[j]*fac;
+    }
+  }
+  row++;
+}
+void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
+                                       const double fac, const Structure &st, const StDescriptors &st_d)
+{
+
+  for (size_t a=0; a<st.natoms(); ++a) {
+    for (size_t jj=0; jj<st_d.fd[a].size(); ++jj) {
+      const size_t j=st.near_neigh_idx[a][jj];
+      const size_t aa = st.get_nn_iindex(a,j,jj);
+      const fd_type &fij = st_d.fd[a][jj];
+      const fd_type &fji = st_d.fd[j][aa];
+      for (size_t k=0; k<3; ++k) {
+        for (size_t d=0; d<fij.rows(); ++d) {
+          Phi(row+k,d) -= fac*(fij(d,k)-fji(d,k));
+        }
+      }
+    }
+    row+=3;
+  }
+}
+void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
+                                        const double fac[6], const Structure &st, const StDescriptors &st_d)
+{
+  for (size_t i=0; i<st.natoms(); ++i) {
+    const Vec3d &ri = st(i).position;
+    for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
+      const size_t j=st.near_neigh_idx[i][jj];
+      const size_t ii = st.get_nn_iindex(i,j,jj);
+      const fd_type &fdij = st_d.fd[i][jj];
+      const fd_type &fdji = st_d.fd[j][ii];
+      const Vec3d &rj = st.nn_pos(i,jj);
+      size_t mn=0;
+      for (size_t x=0; x<3; ++x) {
+        for (size_t y=x; y<3; ++y) {
+          for (size_t d=0; d<fdij.rows(); ++d) {
+            Phi(row+mn,d) += (fdij(d,y)-fdji(d,y))*0.5*fac[mn]*(ri(x)-rj(x));
+          }
+          mn++;
+        }
+      }
+    }
+  }
+  row += 6;
+}

design_matrix/functions/basis_functions/dm_bf_linear.h

+#ifndef DM_BF_LINEAR_H
+#define DM_BF_LINEAR_H
+
+#include "dm_bf_base.h"
+#include "../../../../MODELS/functions/basis_functions/bf_linear.h"
+
+struct DM_BF_Linear: public DM_BF_Base, public BF_Linear
+{
+    DM_BF_Linear();
+    DM_BF_Linear(const Config &c);
+    size_t get_phi_cols(const Config &config);
+    void calc_phi_energy_row(phi_type &Phi, size_t &row,
+            const double fac, const Structure &st, const StDescriptors &st_d);
+    void calc_phi_force_rows(phi_type &Phi, size_t &row,
+            const double fac, const Structure &st, const StDescriptors &st_d);
+    void calc_phi_stress_rows(phi_type &Phi, size_t &row,
+            const double fac[6], const Structure &st, const StDescriptors &st_d);
+};
+#endif

design_matrix/functions/basis_functions/dm_bf_polynomial2.cpp

+#include "dm_bf_polynomial2.h"
+
+Registry<DM_Function_Base>::Register<DM_BF_Polynomial2> DM_BF_Polynomial2_1( "BF_Polynomial2" );
+Registry<DM_Function_Base,Config&>::Register<DM_BF_Polynomial2> DM_BF_Polynomial2_2( "BF_Polynomial2" );
+
+DM_BF_Polynomial2::DM_BF_Polynomial2() {}
+DM_BF_Polynomial2::DM_BF_Polynomial2(const Config &c): BF_Polynomial2(c)
+{}
+size_t DM_BF_Polynomial2::get_phi_cols(const Config &config)
+{
+  size_t cols = config.get<size_t>("DSIZE");
+  return (cols*cols+cols)/2;
+}
+void DM_BF_Polynomial2::calc_phi_energy_row(phi_type &Phi,
+                                            size_t &row,
+                                            const double fac,
+                                            const Structure &,
+                                            const StDescriptors &st_d)
+{
+  for (size_t a=0; a<st_d.naed();++a) {
+    const aed_type2& aed = st_d.get_aed(a);
+    size_t b=0;
+    for (size_t i=0; i<st_d.dim(); ++i) {
+      for (size_t ii=i; ii<st_d.dim(); ++ii) {
+        Phi(row,b++) += aed(i)*aed(ii)*fac;
+      }
+    }
+  }
+  row++;
+}
+void DM_BF_Polynomial2::calc_phi_force_rows(phi_type &Phi,
+                                            size_t &row,
+                                            const double fac,
+                                            const Structure &st,
+                                            const StDescriptors &st_d)
+{
+  for (size_t a=0; a<st.natoms(); ++a) {
+    const aed_type2& aedi = st_d.get_aed(a);
+    for (size_t jj=0; jj<st_d.fd[a].size(); ++jj) {
+      const size_t j=st.near_neigh_idx[a][jj];
+      size_t aa = st.get_nn_iindex(a,j,jj);
+      const fd_type &fdji = st_d.fd[j][aa];
+      const fd_type &fdij = st_d.fd[a][jj];
+      const aed_type2& aedj = st_d.get_aed(j);
+
+      for (size_t k=0; k<3; ++k) {
+        size_t b=0;
+        for (size_t i=0; i<fdij.rows(); ++i) {
+          for (size_t ii=i; ii<fdij.rows(); ++ii) {
+            Phi(row+k,b) -= fac*(fdij(i,k)*aedi(ii) + fdij(ii,k)*aedi(i)
+              - fdji(i,k)*aedj(ii) - fdji(ii,k)*aedj(i));
+            b++;
+          }
+        }
+      }
+
+    }
+    row+=3;
+  }
+}
+void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
+                                             size_t &row,
+                                             const double fac[6],
+                                             const Structure &st, 
+                                             const StDescriptors &st_d)
+{
+  for (size_t a=0; a<st.natoms(); ++a) {
+    const Vec3d &ri = st(a).position;
+    const aed_type2& aedi = st_d.get_aed(a);
+    for (size_t jj=0; jj<st_d.fd[a].size(); ++jj) {
+      const size_t j=st.near_neigh_idx[a][jj];
+      size_t aa = st.get_nn_iindex(a,j,jj);
+      const fd_type &fdji = st_d.fd[j][aa];
+      const fd_type &fdij = st_d.fd[a][jj];
+      const Vec3d &rj = st.nn_pos(a,jj);
+      const aed_type2& aedj = st_d.get_aed(j);
+      size_t mn=0;
+      for (size_t x=0; x<3; ++x) {
+        for (size_t y=x; y<3; ++y) {
+          size_t b=0;
+          for (size_t i=0; i<fdij.rows(); ++i) {
+            for (size_t ii=i; ii<fdij.rows(); ++ii) {
+              Phi(row+mn,b) += 0.5*fac[mn]*(ri(x)-rj(x))
+                *(fdij(i,y)*aedi(ii) + fdij(ii,y)*aedi(i)
+                - fdji(i,y)*aedj(ii) - fdji(ii,y)*aedj(i));
+              b++;
+            }
+          }
+          mn++;
+        }
+      }
+    }
+  }
+  row += 6;
+}

design_matrix/functions/basis_functions/dm_bf_polynomial2.h

+#ifndef DM_BF_POLYNOMIAL2_H
+#define DM_BF_POLYNOMIAL2_H
+#include "dm_bf_base.h"
+#include "../../../../MODELS/functions/basis_functions/bf_polynomial2.h"
+
+struct DM_BF_Polynomial2: public DM_BF_Base, public BF_Polynomial2
+{
+    DM_BF_Polynomial2();
+    DM_BF_Polynomial2(const Config &c);
+    size_t get_phi_cols(const Config &config);
+    void  calc_phi_energy_row(phi_type &Phi, size_t &row,
+            const double fac, const Structure &st, const StDescriptors &st_d);
+    void  calc_phi_force_rows(phi_type &Phi, size_t &row,
+            const double fac, const Structure &st, const StDescriptors &st_d);
+    void  calc_phi_stress_rows(phi_type &Phi, size_t &row,
+            const double fac[6], const Structure &st, const StDescriptors &st_d);
+};
+#endif

design_matrix/functions/dm_function_base.cpp

+#include "dm_function_base.h"
+

design_matrix/functions/dm_function_base.h

+#ifndef DM_FUNCTION_H
+#define DM_FUNCTION_H
+
+#include "../../../CORE/config/config.h"
+#include "../../../CORE/registry.h"
+#include "../../../CORE/core_types.h"
+#include "../../../MODELS/functions/function_base.h"
+#include "../../structure.h"
+#include "../../st_descriptors.h"
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+/** Base class for Kernels and Basis Functions */
+struct DM_Function_Base: public virtual Function_Base {
+
+    // Derived classes must implement Derived() and Derived(Config)
+    virtual ~DM_Function_Base() {};
+
+    virtual size_t get_phi_cols(const Config &)=0;
+    virtual void calc_phi_energy_row(phi_type &, size_t &,
+            const double , const Structure &, const StDescriptors &)=0;
+    virtual void calc_phi_force_rows(phi_type &, size_t &,
+            const double , const Structure &, const StDescriptors &)=0;
+    virtual void calc_phi_stress_rows(phi_type &, size_t &,
+            const double[6], const Structure &, const StDescriptors &)=0;
+};
+template<> inline Registry<DM_Function_Base>::Map Registry<DM_Function_Base>::registry{};
+template<> inline Registry<DM_Function_Base,Config&>::Map Registry<DM_Function_Base,Config&>::registry{};
+#endif