diff --git a/cutoffs/cut_all.h b/cutoffs/cut_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..476c1f830e459c392ff55db9482cb0402855842b
--- /dev/null
+++ b/cutoffs/cut_all.h
@@ -0,0 +1 @@
+#include "cutoffs.h"
diff --git a/cutoffs/cutoffs.cpp b/cutoffs/cutoffs.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3fbd2613fdf13af5537bb069aac374accf53b155
--- /dev/null
+++ b/cutoffs/cutoffs.cpp
@@ -0,0 +1,166 @@
+#include "cutoffs.h"
+#include <cmath>
+#include <stdexcept>
+
+
+Registry<Cut_Base,double>::Register<Cut_Cos> Cut_Cos_1( "Cut_Cos" );
+Registry<Cut_Base,double>::Register<Cut_Tanh> Cut_Tanh_1( "Cut_Tanh" );
+Registry<Cut_Base,double>::Register<Cut_Poly2> Cut_Poly_1( "Cut_Poly2" );
+Registry<Cut_Base,double>::Register<Cut_Dummy> Cut_Dummy_1( "Cut_Dummy" );
+
+
+Cut_Base::~Cut_Base() {}
+void Cut_Base::test_rcut(const double r)
+{
+ if (r<0.0)
+ throw std::runtime_error("Cutoff distance cannot be negative.");
+}
+
+Cut_Cos::Cut_Cos() {}
+Cut_Cos::Cut_Cos(double rcut)
+{
+ set_rcut(rcut);
+ test_rcut(rcut);
+}
+
+std::string Cut_Cos::label() {
+ return lab;
+}
+
+void Cut_Cos::set_rcut(const double r) {
+ test_rcut(r);
+ rcut=r;
+ rcut_sq=r*r;
+ rcut_inv= r<=0 ? 0.0 : 1.0/r;
+}
+double Cut_Cos::get_rcut() {
+ return rcut;
+}
+
+double Cut_Cos::get_rcut_sq() {
+ return rcut_sq;
+}
+
+double Cut_Cos::calc(double r) {
+ if (r>=rcut) return 0.0;
+ return 0.5 * (std::cos(M_PI*r*rcut_inv) + 1.0);
+}
+double Cut_Cos::calc_prime(double r) {
+ if (r>=rcut) return 0.0;
+ double t = M_PI*rcut_inv;
+ return -0.5 * t*sin(t*r);
+
+}
+
+Cut_Dummy::Cut_Dummy() {}
+Cut_Dummy::Cut_Dummy(double rcut)
+{
+ set_rcut(rcut),
+ test_rcut(rcut);
+}
+
+std::string Cut_Dummy::label() {
+ return lab;
+}
+
+void Cut_Dummy::set_rcut(const double r) {
+ test_rcut(r);
+ rcut=r;
+ rcut_sq=r*r;
+ rcut_inv= r<=0 ? 0.0 : 1.0/r;
+}
+double Cut_Dummy::get_rcut() {
+ return rcut;
+}
+
+double Cut_Dummy::get_rcut_sq() {
+ return rcut_sq;
+}
+
+double Cut_Dummy::calc(double r) {
+ if (r>=rcut) return 0.0;
+ return 1.0;
+}
+double Cut_Dummy::calc_prime(double) {
+ return 0.0;
+}
+
+Cut_Tanh::Cut_Tanh() {}
+Cut_Tanh::Cut_Tanh(double rcut)
+{
+ set_rcut(rcut);
+ test_rcut(rcut);
+}
+
+std::string Cut_Tanh::label() {
+ return lab;
+}
+
+void Cut_Tanh::set_rcut(const double r) {
+ test_rcut(r);
+ rcut=r;
+ rcut_sq=r*r;
+ rcut_inv= r<=0 ? 0.0 : 1.0/r;
+}
+double Cut_Tanh::get_rcut() {
+ return rcut;
+}
+
+double Cut_Tanh::get_rcut_sq() {
+ return rcut_sq;
+}
+
+double Cut_Tanh::calc(double r) {
+ if (r>=rcut) return 0.0;
+ double t = std::tanh(1.0 - r*rcut_inv);
+ return t*t*t;
+}
+double Cut_Tanh::calc_prime(double r) {
+ if (r>=rcut) return 0.0;
+ double a = (1.0-r*rcut_inv);
+ double t = std::tanh(a);
+ double s = 1.0/(std::cosh(a));
+ return -3.0*t*t*rcut_inv*s*s;
+}
+
+Cut_Poly2::Cut_Poly2() {}
+Cut_Poly2::Cut_Poly2(double rcut)
+{
+ if (rcut<=1.0)
+ throw std::runtime_error(
+ "Cut_Poly2 cutoff distance cannot be smaller than 1.0.");
+ set_rcut(rcut);
+ test_rcut(rcut);
+}
+
+std::string Cut_Poly2::label() {
+ return lab;
+}
+
+void Cut_Poly2::set_rcut(const double r) {
+ test_rcut(r);
+ rcut=r;
+ rcut_sq=r*r;
+ rcut_inv= r<=0 ? 0.0 : 1.0/r;
+ rcut_inner=rcut-1.0;
+}
+double Cut_Poly2::get_rcut() {
+ return rcut;
+}
+
+double Cut_Poly2::get_rcut_sq() {
+ return rcut_sq;
+}
+
+double Cut_Poly2::calc(double r) {
+ if (r>=rcut) return 0.0;
+ else if (r<= rcut_inner) return 1.0;
+ double rs=r-rcut_inner;
+ return rs*rs*rs*(rs*(15.0-6.0*rs)-10.0)+1;
+}
+double Cut_Poly2::calc_prime(double r) {
+ if (r>=rcut) return 0.0;
+ else if (r<= rcut_inner) return 0.0;
+ double rs=r-rcut_inner;
+ return -30.0*(rs-1.0)*(rs-1.0)*rs*rs;
+}
diff --git a/cutoffs/cutoffs.h b/cutoffs/cutoffs.h
new file mode 100644
index 0000000000000000000000000000000000000000..98c5622dd1374e6348c2cecdc29ace547c36a2c0
--- /dev/null
+++ b/cutoffs/cutoffs.h
@@ -0,0 +1,143 @@
+#ifndef CUTOFFS_H
+#define CUTOFFS_H
+
+#include <string>
+#include "../../CORE/config/config.h"
+#include "../../CORE/registry.h"
+
+/** Base class to be inherited by all cutoffs */
+class Cut_Base {
+ public:
+ virtual ~Cut_Base();
+ virtual std::string label()=0;
+ virtual double calc(double r)=0;
+ virtual double get_rcut()=0;
+ virtual void set_rcut(const double r)=0;
+ virtual double get_rcut_sq()=0;
+ virtual double calc_prime(double r)=0;
+ void test_rcut(const double r);
+};
+
+/** \brief Dummy cutoff function
+ * \f[
+ * f_c(r) =
+ * \begin{equation}
+ * \begin{cases}
+ * 1 & \text{if } r \leq r_c\\
+ * 0 & \text{otherwise}\\
+ * \end{cases}
+ * \end{equation}
+ * \f]
+ */
+class Cut_Dummy : public Cut_Base {
+ private:
+ std::string lab = "Cut_Dummy";
+ double rcut, rcut_sq, rcut_inv;
+ public:
+ Cut_Dummy();
+ Cut_Dummy(double rcut);
+ std::string label() ;
+ void set_rcut(const double r);
+ double get_rcut();
+ double get_rcut_sq();
+ double calc(double r);
+ double calc_prime(double r);
+};
+
+/** \brief Cos cutoff function
+ *
+ * \f[
+ * f_c(r) =
+ * \begin{equation}
+ * \begin{cases}
+ * \frac{1}{2}\big[ \cos\big( \frac{\pi r}{r_c} \big)+1 \big] & \text{if } r \leq r_c\\
+ * 0 & \text{otherwise}\\
+ * \end{cases}
+ * \end{equation}
+ * \f]
+ *
+ * <div class="csl-entry">Behler, J., Parrinello, M. (2007).
+ * Generalized neural-network representation of high-dimensional
+ * potential-energy surfaces. <i>Physical Review Letters</i>,
+ * <i>98</i>(14), 146401. https://doi.org/10.1103/PhysRevLett.98.146401</div>
+ */
+class Cut_Cos : public Cut_Base {
+ private:
+ std::string lab = "Cut_Cos";
+ double rcut, rcut_sq, rcut_inv;
+ public:
+ Cut_Cos();
+ Cut_Cos(double rcut);
+ std::string label() ;
+ void set_rcut(const double r);
+ double get_rcut();
+ double get_rcut_sq();
+ double calc(double r);
+ double calc_prime(double r);
+};
+
+/** \brief Tanh cutoff function.
+ * \f[
+ * f_c(r) =
+ * \begin{equation}
+ * \begin{cases}
+ * \tanh^3\big( 1 -\frac{r}{r_c} \big) & \text{if } r \leq r_c\\
+ * 0 & \text{otherwise}\\
+ * \end{cases}
+ * \end{equation}
+ * \f]
+ *
+ * <div class="csl-entry">Behler, J. (2011). Atom-centered symmetry
+ * functions for constructing high-dimensional neural network potentials.
+ * <i>J. Chem. Phys.</i>, <i>134</i>(7), 074106.
+ * https://doi.org/10.1063/1.3553717</div>
+ */
+class Cut_Tanh : public Cut_Base {
+ private:
+ std::string lab = "Cut_Tanh";
+ double rcut, rcut_sq, rcut_inv;
+ public:
+ Cut_Tanh();
+ Cut_Tanh(double rcut);
+ std::string label() ;
+ void set_rcut(const double r);
+ double get_rcut();
+ double get_rcut_sq();
+ double calc(double r);
+ double calc_prime(double r);
+};
+
+/** \brief Polynomial-2 cutoff function.
+ *
+ * \f[
+ * f_c(r) =
+ * \begin{equation}
+ * \begin{cases}
+ * 1 & \text{if } r \leq (r_c-1)\\
+ * r^3(r(15-6r)-10)+1 & \text{if } (r_c-1) < r \leq r_c\\
+ * 0 & \text{otherwise}\\
+ * \end{cases}
+ * \end{equation}
+ * \f]
+
+ *<div class="csl-entry">Singraber, A., Rg Behler, J., Dellago, C. (2019).
+ Library-Based LAMMPS Implementation of High-Dimensional
+ Neural Network Potentials. https://doi.org/10.1021/acs.jctc.8b00770</div>
+ */
+class Cut_Poly2 : public Cut_Base {
+ private:
+ std::string lab = "Cut_Poly2";
+ double rcut, rcut_sq, rcut_inv;
+ double rcut_inner;
+ public:
+ Cut_Poly2();
+ Cut_Poly2(double rcut);
+ std::string label() ;
+ void set_rcut(const double r);
+ double get_rcut();
+ double get_rcut_sq();
+ double calc(double r);
+ double calc_prime(double r);
+};
+template<> inline Registry<Cut_Base,double>::Map Registry<Cut_Base,double>::registry{};
+#endif
diff --git a/dc_selector.h b/dc_selector.h
new file mode 100644
index 0000000000000000000000000000000000000000..a51cf8e3348059975b00911fcdc01fbdd17b6c3a
--- /dev/null
+++ b/dc_selector.h
@@ -0,0 +1,116 @@
+#ifndef DC_SELECTOR_H
+#define DC_SELECTOR_H
+
+#include "../CORE/config/config.h"
+#include "descriptors/d_all.h"
+#include "cutoffs/cut_all.h"
+#include "../CORE/registry.h"
+
+class DC_Selector {
+ public:
+ Cut_Base *c2b=nullptr;
+ Cut_Base *c3b=nullptr;
+ Cut_Base *cmb=nullptr;
+ D2_Base *d2b=nullptr;
+ D3_Base *d3b=nullptr;
+ DM_Base *dmb=nullptr;
+
+ Config config;
+
+ DC_Selector ()
+ {};
+
+ // Copy constructor
+ DC_Selector(const DC_Selector& )
+ {
+ std::cout << "Copy constructor called\n";
+ }
+ // Const Assignment operator
+ // Instead of deep copying
+ // we use a trick where we just copy
+ // the config file and run init()
+ // as in a costructor
+ DC_Selector& operator=(const DC_Selector& dc)
+ {
+ config=dc.config;
+ init();
+ return *this;
+ }
+ // Assignment operator
+ DC_Selector& operator=(DC_Selector& dc)
+ {
+ std::swap(config,dc.config);
+ std::swap(c2b,dc.c2b);
+ std::swap(c3b,dc.c3b);
+ std::swap(cmb,dc.cmb);
+ std::swap(d2b,dc.d2b);
+ std::swap(d3b,dc.d3b);
+ std::swap(dmb,dc.dmb);
+ return *this;
+ }
+
+ DC_Selector (const Config &c):
+ config(c)
+ {
+ init();
+ };
+ void init() {
+ double rcutzero = 0.0; // for dummies
+
+ size_t bias=0;
+ if (config.get<bool>("BIAS"))
+ bias++;
+
+ if (config.get<bool>("INIT2B")) {
+ double rcut2b = config.get<double>("RCUT2B");
+ c2b = factory<Cut_Base,double>( config.get<std::string>("RCTYPE2B"), rcut2b );
+ d2b = factory<D2_Base,Config&>( config.get<std::string>("TYPE2B"), config );
+ d2b->fidx = bias;
+ }
+ else {
+ c2b = factory<Cut_Base,double>( "Cut_Dummy", rcutzero );
+ d2b = factory<D2_Base,Config&>( "D2_Dummy", config );
+ d2b->fidx = bias;
+ }
+
+ if (config.get<bool>("INIT3B")) {
+ double rcut3b = config.get<double>("RCUT3B");
+ c3b = factory<Cut_Base,double>( config.get<std::string>("RCTYPE3B"), rcut3b );
+ d3b = factory<D3_Base,Config&>( config.get<std::string>("TYPE3B"), config );
+ // d3b->fidx = bias;
+ }
+ else {
+ c3b = factory<Cut_Base,double>( "Cut_Dummy", rcutzero );
+ d3b = factory<D3_Base,Config&>( "D3_Dummy", config );
+ // d3b->fidx = bias;
+ }
+
+ if (config.get<bool>("INITMB")) {
+ double rcutmb = config.get<double>("RCUTMB");
+ cmb = factory<Cut_Base,double>( config.get<std::string>("RCTYPEMB"), rcutmb );
+ dmb = factory<DM_Base,Config&>( config.get<std::string>("TYPEMB"), config );
+ dmb->fidx = bias + d2b->size();
+ }
+ else {
+ cmb = factory<Cut_Base,double>( "Cut_Dummy", rcutzero );
+ dmb = factory<DM_Base,Config&>( "DM_Dummy", config );
+ dmb->fidx = bias + d2b->size();
+ }
+ }
+ ~DC_Selector()
+ {
+ if (d2b)
+ delete d2b;
+ if (c2b)
+ delete c2b;
+ if (c3b)
+ delete c3b;
+ if (d3b)
+ delete d3b;
+ if (cmb)
+ delete cmb;
+ if (dmb)
+ delete dmb;
+ }
+};
+#endif
diff --git a/descriptors/d2/d2_all.h b/descriptors/d2/d2_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..684650372ef5ea238567731784c9a4aa317a386a
--- /dev/null
+++ b/descriptors/d2/d2_all.h
@@ -0,0 +1,6 @@
+#include "d2_base.h"
+#include "d2_dummy.h"
+#include "d2_blip.h"
+#include "d2_bp.h"
+#include "d2_eam.h"
+#include "d2_lj.h"
diff --git a/descriptors/d2/d2_base.cpp b/descriptors/d2/d2_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..017b0ce5223745e7e82d95a1066e7fdb14fe2650
--- /dev/null
+++ b/descriptors/d2/d2_base.cpp
@@ -0,0 +1,14 @@
+#include "d2_base.h"
+//template struct Registry<D2_Base, Config &>;
+void D2_Base::calc_fd_approx(
+ const double ,
+ const double ,
+ const double ,
+ aed_rtype aedp,
+ aed_rtype aedn,
+ fd_type &fd_ij,
+ const double h) {
+
+ for (size_t i=0; i<(*this).size(); ++i)
+ fd_ij(i,0) = (aedp(i)-aedn(i))/(2*h);
+}
diff --git a/descriptors/d2/d2_base.h b/descriptors/d2/d2_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..6ff85dba928f04fd656c1ba2c5b440eb1d3a6784
--- /dev/null
+++ b/descriptors/d2/d2_base.h
@@ -0,0 +1,68 @@
+#ifndef D2_BASE_H
+#define D2_BASE_H
+
+#include "../d_base.h"
+
+/** \brief Base class for all two-body type descriptors.
+ *
+ * All two-body descriptors **must** inherit this class.
+ */
+class D2_Base: public D_Base {
+ public:
+ size_t fidx=0; // first index in aed and fd arrays for this descriptor
+ virtual ~D2_Base() {};
+
+ /** \brief Calculate \ref AED
+ *
+ * Calculate Atomic Energy Descriptor for the atom local environment.
+ */
+ virtual void calc_aed(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ aed_rtype aed)=0;
+
+ /** \brief Calculate \ref FD
+ *
+ * Calculate Force Descriptor for the atom local environment.
+ */
+ virtual void calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij)=0;
+
+ /** \brief Calculate approximation to \ref FD
+ *
+ * Calculate Force Descriptor for the atom local environment
+ * using central difference approximation.
+ *
+ * h -> something small e.g. 1e-8
+ * fcijp -> cutoff(rij+h,...)
+ * fcijn -> cutoff(rij-h,...)
+ * aedp -> calc_aed(rij+h,...)
+ * aedn -> calc_aed(rij-h,...)
+ * fd_ij - descriptor derivative
+ */
+ virtual void calc_fd_approx(
+ const double rij,
+ const double fcijp,
+ const double fcijn,
+ aed_rtype aedp,
+ aed_rtype aedn,
+ fd_type &fd_ij,
+ const double h);
+
+ /** \brief Calculate \ref AED + \ref FD */
+ virtual void calc_all(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij)=0;
+ virtual std::string label()=0;
+};
+template<> inline Registry<D2_Base,Config&>::Map Registry<D2_Base,Config&>::registry{};
+#endif
diff --git a/descriptors/d2/d2_blip.cpp b/descriptors/d2/d2_blip.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..35640daf64e59f191ca18fb904691791147ef75d
--- /dev/null
+++ b/descriptors/d2/d2_blip.cpp
@@ -0,0 +1,82 @@
+#include "d2_blip.h"
+#include "../d_basis_functions.h"
+#include <stdexcept>
+
+Registry<D2_Base,Config&>::Register<D2_Blip> D2_Blip_1( "D2_Blip" );
+
+D2_Blip::D2_Blip(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{
+
+ if (!c.get<bool>("INIT2B")) return;
+
+ get_grid(c,"CGRID2B",mius);
+ get_grid(c,"SGRID2B",etas);
+
+ if (verbose) {
+ std::cout << std::endl;
+ std::cout << "SGRID2B: " << etas.size() << std::endl;
+ for (auto e:etas) std::cout << e << " ";
+
+ std::cout << std::endl;
+ std::cout << "CGRID2B: " << mius.size() << std::endl;
+ for (auto m:mius) std::cout << m << " ";
+ std::cout << std::endl;
+ }
+
+ if (mius.size()!=etas.size()) {
+ throw std::runtime_error("SGRID2B and CGRID2B arrays differ in size.\n");
+ }
+
+ s=mius.size();
+}
+
+void D2_Blip::calc_aed(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ aed_rtype aed)
+{
+
+ size_t i=fidx;
+ for (size_t c=0; c<mius.size(); c++) {
+ double t = B(etas[c]*(rij-mius[c]),fc_ij);
+ aed(i++) += t;
+ }
+
+}
+void D2_Blip::calc_dXijdri(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij)
+{
+ size_t i=fidx;
+ for (size_t c=0; c<mius.size(); c++) {
+ fd_ij(i++,0) = dB(etas[c]*(rij-mius[c]),etas[c],fc_ij,fcp_ij);
+ }
+}
+void D2_Blip::calc_all(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij)
+{
+
+ size_t i=fidx;
+ for (size_t c=0; c<mius.size(); c++) {
+ aed(i) += B(etas[c]*(rij-mius[c]),fc_ij);
+ fd_ij(i,0) = dB(etas[c]*(rij-mius[c]),etas[c],fc_ij,fcp_ij);
+ ++i;
+ }
+}
+
+size_t D2_Blip::size() {
+ return s;
+}
+std::string D2_Blip::label() {
+ return lab;
+}
diff --git a/descriptors/d2/d2_blip.h b/descriptors/d2/d2_blip.h
new file mode 100644
index 0000000000000000000000000000000000000000..beb7d762af4d58c3924fb8932cc4443e3eb763be
--- /dev/null
+++ b/descriptors/d2/d2_blip.h
@@ -0,0 +1,77 @@
+#ifndef D2_BLIP_H
+#define D2_BLIP_H
+#include "d2_base.h"
+#include <vector>
+
+/** \brief Blip two-body descriptor.
+ *
+ * \f[
+ * V_i^{\eta,r_s} =\sum_{j \neq i} B(\eta(r_{ij}-r_s))f_c(r_{ij})
+ * \f]
+ *
+ * where \f$ f_c \f$ is a cutoff function and \f$ B \f$ is a blip
+ * basis function centered at \f$r_s\f$ of width \f$4/\eta\f$.
+ *
+ * \ref CGRID2B parameters control position \f$ r_s \f$ of blip centres.
+ *
+ * \ref SGRID2B parameters control width \f$ \eta \f$ of blips.
+ *
+ * Blip basis function is built out of 3rd degree polynomials
+ * in the four intervals [-2,-1], [-1,0], [0,1], [1,2] and is defined as:
+ *
+ * \f[
+ \begin{equation}
+ B(r) =
+ \begin{cases}
+ 1-\frac{3}{2}r^2+\frac{3}{4}|r|^3 & \text{if} \qquad 0<|r|<1\\
+ \frac{1}{4}(2-|r|)^3 & \text{if} \qquad 1<|r|<2\\
+ 0 & \text{if} \qquad |r|>2
+ \end{cases}
+ \end{equation}
+ * \f]
+ *
+ * More details about the blip basis functions can be found in the following paper:
+ *
+ * <div class="csl-entry">Hernández, E., Gillan, M., Goringe, C.
+ * (1997). Basis functions for linear-scaling first-principles calculations.
+ * <i>Physical Review B - Condensed Matter and Materials Physics</i>,
+ * <i>55</i>(20), 13485–13493. https://doi.org/10.1103/PhysRevB.55.13485</div>
+ *
+ * Required keys:
+ * \ref INIT2B \ref CGRID2B \ref SGRID2B
+ */
+class D2_Blip : public D2_Base {
+ private:
+ size_t s=0;
+ std::string lab="D2_Blip";
+ v_type etas;
+ v_type mius;
+ double get_blip(double);
+ double get_dblip(double, double);
+ v_type gen_blip_grid(double, double);
+ int verbose;
+
+ public:
+ D2_Blip(Config &config);
+ void calc_aed(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ aed_rtype aed);
+ void calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij);
+ void calc_all(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d2/d2_bp.cpp b/descriptors/d2/d2_bp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..abdce9df95010d09363ad5bfb2b434a8d713745c
--- /dev/null
+++ b/descriptors/d2/d2_bp.cpp
@@ -0,0 +1,79 @@
+#include "d2_bp.h"
+#include "../d_basis_functions.h"
+
+Registry<D2_Base,Config&>::Register<D2_BP> D2_BP_1( "D2_BP" );
+
+D2_BP::D2_BP(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{
+
+ if (!c.get<bool>("INIT2B")) return;
+
+ get_grid(c,"CGRID2B",mius);
+ get_grid(c,"SGRID2B",etas);
+
+ if (verbose) {
+ std::cout << std::endl;
+ std::cout << "SGRID2B: " << etas.size() << std::endl;
+ for (auto e:etas) std::cout << e << " ";
+ std::cout << std::endl;
+
+ std::cout << "CGRID2B: " << mius.size() << std::endl;
+ for (auto m:mius) std::cout << m << " ";
+ std::cout << std::endl;
+ }
+
+ if (mius.size()!=etas.size()) {
+ throw std::runtime_error("SGRID2B and CGRID2B arrays differ in size.\n");
+ }
+
+ s=mius.size();
+}
+
+void D2_BP::calc_aed(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ aed_rtype aed)
+{
+
+ size_t i=fidx;
+ for (size_t c=0; c<mius.size(); c++) {
+ aed(i++) += G(rij,etas[c],mius[c],fc_ij);
+ }
+
+}
+void D2_BP::calc_dXijdri(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij)
+{
+ size_t i=fidx;
+ for (size_t c=0; c<mius.size(); c++) {
+ fd_ij(i++,0) = dG(rij,etas[c],mius[c],fc_ij,fcp_ij);
+ }
+}
+void D2_BP::calc_all(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij)
+{
+
+ size_t i=fidx;
+ for (size_t c=0; c<mius.size(); c++) {
+ aed(i) += G(rij,etas[c],mius[c],fc_ij);
+ fd_ij(i,0) = dG(rij,etas[c],mius[c],fc_ij,fcp_ij);
+ ++i;
+ }
+}
+size_t D2_BP::size() {
+ return s;
+}
+std::string D2_BP::label() {
+ return lab;
+}
diff --git a/descriptors/d2/d2_bp.h b/descriptors/d2/d2_bp.h
new file mode 100644
index 0000000000000000000000000000000000000000..d73e87371d518ef8747a97ba056ecbb5757d0a3e
--- /dev/null
+++ b/descriptors/d2/d2_bp.h
@@ -0,0 +1,58 @@
+#ifndef D2_BP_H
+#define D2_BP_H
+#include "d2_base.h"
+#include <vector>
+
+/** \brief Behler-Parrinello two-body descriptor.
+ *
+ * \f[
+ * V_i^{\eta,r_s} = \sum_{j \neq i} \exp{\Big(-\eta(r_{ij}-r_s)^2\Big)}f_c(r_{ij})
+ * \f]
+ *
+ * \ref CGRID2B parameters control position \f$ r_s \f$ of the gaussian basis function.
+ *
+ * \ref SGRID2B parameters control width \f$ \eta \f$ of the gaussian basis function.
+ *
+ * This is essentially a \f$ G^1_i \f$ descriptor from the below paper with
+ * an exception that it can use any cutoff function defined in Ta-dah!:
+ *
+ * <div class="csl-entry">Behler, J., Parrinello, M. (2007).
+ * Generalized neural-network representation of high-dimensional
+ * potential-energy surfaces. <i>Physical Review Letters</i>,
+ * <i>98</i>(14), 146401. https://doi.org/10.1103/PhysRevLett.98.146401</div>
+ *
+ * Required Config keys:
+ * \ref INIT2B \ref CGRID2B \ref SGRID2B
+ */
+class D2_BP : public D2_Base {
+ private:
+ size_t s=0;
+ std::string lab="D2_BP";
+ v_type etas;
+ v_type mius;
+ int verbose;
+
+ public:
+ D2_BP(Config &config);
+ void calc_aed(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ aed_rtype aed);
+ void calc_dXijdri(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij);
+ void calc_all(
+ const double rij,
+ const double ,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d2/d2_dummy.cpp b/descriptors/d2/d2_dummy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ea73bac21fe1eaef4706da13aee7eb353aa14b9
--- /dev/null
+++ b/descriptors/d2/d2_dummy.cpp
@@ -0,0 +1,29 @@
+#include "d2_dummy.h"
+
+Registry<D2_Base,Config&>::Register<D2_Dummy> D2_Dummy_1( "D2_Dummy" );
+
+D2_Dummy::D2_Dummy() {}
+D2_Dummy::D2_Dummy(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{}
+
+void D2_Dummy::calc_aed(
+ const double,
+ const double,
+ const double,
+ aed_rtype ) {}
+void D2_Dummy::calc_dXijdri(
+ const double,
+ const double,
+ const double,
+ const double,
+ fd_type &) {}
+void D2_Dummy::calc_all(
+ const double,
+ const double,
+ const double,
+ const double,
+ aed_rtype ,
+ fd_type &) {}
+size_t D2_Dummy::size() { return s; }
+std::string D2_Dummy::label() { return lab; }
diff --git a/descriptors/d2/d2_dummy.h b/descriptors/d2/d2_dummy.h
new file mode 100644
index 0000000000000000000000000000000000000000..87a73fd640747397e04ce4ea42c1d1194c3a4b0a
--- /dev/null
+++ b/descriptors/d2/d2_dummy.h
@@ -0,0 +1,41 @@
+#ifndef D2_DUMMY_H
+#define D2_DUMMY_H
+#include "d2_base.h"
+
+/**
+ * Dummy two-body descriptor.
+ *
+ * Use it to satisfy DescriptorsCalc requirements in case
+ * when two-body descriptor is not required.
+ *
+ */
+class D2_Dummy : public D2_Base {
+ private:
+ size_t s=0;
+ std::string lab="D2_Dummy";
+ int verbose;
+ public:
+ D2_Dummy();
+ D2_Dummy(Config &);
+ void calc_aed(
+ const double,
+ const double,
+ const double,
+ aed_rtype aed);
+ void calc_dXijdri(
+ const double,
+ const double,
+ const double,
+ const double,
+ fd_type &fd_ij);
+ void calc_all(
+ const double,
+ const double,
+ const double,
+ const double,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d2/d2_eam.cpp b/descriptors/d2/d2_eam.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..99bcbf14837cb16641cc83133e0c88dfabe989cb
--- /dev/null
+++ b/descriptors/d2/d2_eam.cpp
@@ -0,0 +1,176 @@
+#include "d2_eam.h"
+
+Registry<D2_Base,Config&>::Register<D2_EAM> D2_EAM_1( "D2_EAM" );
+
+D2_EAM::D2_EAM(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{
+ if (!c.get<bool>("INIT2B")) {
+ s=0;
+ return;
+ }
+ ef.file_path = c("SETFL")[0];
+ read_setfl();
+
+ if (abs(ef.rcut - c.get<double>("RCUT2B")) > std::numeric_limits<double>::min()) {
+ if (verbose) {
+ std::cout << std::endl;
+ std::cout << "Config file cutoff and setfl file cutoff differ: "
+ << c.get<double>("RCUT2B") << " " << ef.rcut << std::endl;
+ std::cout << "Enforcing SETFL file cutoff: " << ef.rcut << std::endl;
+ }
+ c.remove("RCUT2B");
+ c.add("RCUT2B", ef.rcut);
+ c.postprocess();
+ }
+
+
+ frho_spline.resize(ef.nrho+1, std::vector<double>(7));
+ rhor_spline.resize(ef.nr+1, std::vector<double>(7));
+ z2r_spline.resize(ef.nr+1, std::vector<double>(7));
+
+
+ gen_splines(ef.nrho, ef.drho, ef.frho, frho_spline);
+ gen_splines(ef.nr, ef.dr, ef.rhor, rhor_spline);
+ gen_splines(ef.nr, ef.dr, ef.z2r, z2r_spline);
+
+}
+
+void D2_EAM::calc_aed(
+ const double rij,
+ const double,
+ const double,
+ aed_rtype aed)
+{
+
+ double r = rij;
+ const double recip = 1.0/r;
+ double p = r*ef.rdr + 1.0;
+ int m = static_cast<int> (p);
+ m = std::min(m,ef.nr-1);
+ p -= m;
+ p = std::min(p,1.0);
+ double z2 = ((z2r_spline[m][3]*p + z2r_spline[m][4])*p + z2r_spline[m][5])*p + z2r_spline[m][6];
+ aed(fidx) += z2*recip;
+
+}
+void D2_EAM::calc_dXijdri(
+ const double rij,
+ const double,
+ const double,
+ const double,
+ fd_type &fd_ij)
+{
+
+ const double r = rij;
+ const double recip = 1.0/r;
+ double p = r*ef.rdr + 1.0;
+ int m = static_cast<int> (p);
+ m = std::min(m,ef.nr-1);
+ p -= m;
+ p = std::min(p,1.0);
+ double z2p = (z2r_spline[m][0]*p + z2r_spline[m][1])*p + z2r_spline[m][2];
+ double z2 = ((z2r_spline[m][3]*p + z2r_spline[m][4])*p + z2r_spline[m][5])*p + z2r_spline[m][6];
+ //double phi = z2*recip;
+ //double phip = z2p*recip - z2*recip*recip;
+
+ fd_ij(fidx,0) = (z2p*recip - z2*recip*recip);
+ //force_fp_ij[first_idx] = 0.5*(z2p*recip - z2*recip*recip);
+
+}
+void D2_EAM::calc_all(
+ const double rij,
+ const double,
+ const double,
+ const double,
+ aed_rtype aed,
+ fd_type &fd_ij)
+{
+ double r = rij;
+ const double recip = 1.0/r;
+ double p = r*ef.rdr + 1.0;
+ int m = static_cast<int> (p);
+ m = std::min(m,ef.nr-1);
+ p -= m;
+ p = std::min(p,1.0);
+ double z2 = ((z2r_spline[m][3]*p + z2r_spline[m][4])*p + z2r_spline[m][5])*p + z2r_spline[m][6];
+ aed(fidx) += z2*recip;
+
+ double z2p = (z2r_spline[m][0]*p + z2r_spline[m][1])*p + z2r_spline[m][2];
+ // 0.5 b/c full neighbour list is used TODO check
+ fd_ij(fidx,0) = (z2p*recip - z2*recip*recip);
+
+
+}
+size_t D2_EAM::size() {
+ return s;
+}
+std::string D2_EAM::label() {
+ return lab;
+}
+
+void D2_EAM::read_setfl()
+{
+ std::string line;
+ std::ifstream in_file(ef.file_path);
+ if (!in_file.good())
+ throw std::runtime_error("SETFL file does not exists.\n");
+
+ if (in_file.is_open()) {
+ if (verbose)
+ std::cout << std::endl << "<D2_EAM> Reading setfl: " << ef.file_path << std::endl;
+ // skip ficgridt 3 comment lines
+ getline (in_file,line);
+ getline (in_file,line);
+ getline (in_file,line);
+ // skip number of types and types
+ getline (in_file,line);
+ // read 5th line
+ in_file >> ef.nrho >> ef.drho >> ef.nr >> ef.dr >> ef.rcut;
+ in_file >> ef.atomic_number >> ef.atomic_mass >> ef.lattice_param >> ef.lattice;
+ ef.rdr = 1.0/ef.dr;
+ ef.rdrho = 1.0/ef.drho;
+ // prepare arrays
+ ef.frho.resize(ef.nrho);
+ ef.rhor.resize(ef.nr);
+ ef.z2r.resize(ef.nr);
+ // read all data
+ for (int i=0; i<ef.nrho; ++i) in_file >> ef.frho[i];
+ for (int i=0; i<ef.nr; ++i) in_file >> ef.rhor[i];
+ for (int i=0; i<ef.nr; ++i) in_file >> ef.z2r[i];
+ in_file.close();
+ }
+ else {
+ if (verbose) std::cout << "<D2_EAM> Unable to open file: " << ef.file_path << std::endl;;
+ }
+}
+void D2_EAM::gen_splines(int &n, double &delta, std::vector<double> &f,
+ std::vector<std::vector<double>> &spline)
+{
+ // in lammps f is n+1, here is size n
+ for (int m=1; m<=n; m++) spline[m][6] = f[m-1];
+
+ spline[1][5] = spline[2][6] - spline[1][6];
+ spline[2][5] = 0.5 * (spline[3][6]-spline[1][6]);
+ spline[n-1][5] = 0.5 * (spline[n][6]-spline[n-2][6]);
+ spline[n][5] = spline[n][6] - spline[n-1][6];
+
+ for (int m = 3; m <= n-2; m++)
+ spline[m][5] = ((spline[m-2][6]-spline[m+2][6]) +
+ 8.0*(spline[m+1][6]-spline[m-1][6])) / 12.0;
+
+ for (int m = 1; m <= n-1; m++) {
+ spline[m][4] = 3.0*(spline[m+1][6]-spline[m][6]) - 2.0*spline[m][5] - spline[m+1][5];
+ spline[m][3] = spline[m][5] + spline[m+1][5] - 2.0*(spline[m+1][6]-spline[m][6]);
+ }
+
+ spline[n][4] = 0.0;
+ spline[n][3] = 0.0;
+
+ for (int m = 1; m <= n; m++) {
+ spline[m][2] = spline[m][5]/delta;
+ spline[m][1] = 2.0*spline[m][4]/delta;
+ spline[m][0] = 3.0*spline[m][3]/delta;
+ }
+}
+
diff --git a/descriptors/d2/d2_eam.h b/descriptors/d2/d2_eam.h
new file mode 100644
index 0000000000000000000000000000000000000000..6a8545d4b33e1404df48af0e6fc8a37a71cae56b
--- /dev/null
+++ b/descriptors/d2/d2_eam.h
@@ -0,0 +1,79 @@
+#ifndef D2_EAM_H
+#define D2_EAM_H
+#include "d2_base.h"
+
+/** \brief Pair-wise part for the Embedded Atom Method descriptor.
+ *
+ * \f[
+ * V_i = \frac{1}{2} \sum_{j \neq i} \psi(r_{ij})
+ * \f]
+ *
+ * This descriptor will load tabulated values for the two-body potential \f$ \phi \f$
+ * from the provided \ref SETFL file.
+ *
+ * This descriptor is usually used together with the many-body descriptor \ref DM_EAM
+ * although this is not required and user can mix it with any other descriptors
+ * or use it on its own.
+ *
+ * This descriptor will enforce cutoff distance as specified in a \ref SETFL file.
+ * Set \ref RCUT2B to the same value to suppress the warning message.
+ *
+ * Required Config keys:
+ * \ref INIT2B \ref SETFL
+ */
+class D2_EAM : public D2_Base {
+ private:
+ struct eam_file {
+ std::string file_path;
+ std::vector<double> frho;
+ std::vector<double> rhor;
+ std::vector<double> z2r;
+ int nrho=0;
+ double drho=0;
+ int nr;
+ double dr;
+ double rdr;
+ double rdrho;
+ double rcut;
+ int atomic_number;
+ double atomic_mass;
+ double lattice_param;
+ std::string lattice;
+ };
+ eam_file ef;
+
+ std::vector<std::vector<double>> frho_spline;
+ std::vector<std::vector<double>> rhor_spline;
+ std::vector<std::vector<double>> z2r_spline;
+
+ int verbose;
+
+ size_t s=1;
+ std::string lab="D2_EAM";
+ void read_setfl();
+ void gen_splines(int &n, double &delta, std::vector<double> &f, std::vector<std::vector<double>> &spline);
+
+ public:
+ D2_EAM(Config &config);
+ void calc_aed(
+ double rij,
+ const double rij_sq,
+ const double,
+ aed_rtype aed);
+ void calc_dXijdri(
+ double rij,
+ const double rij_sq,
+ const double,
+ const double,
+ fd_type &fd_ij);
+ void calc_all(
+ double rij,
+ const double rij_sq,
+ const double,
+ const double,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d2/d2_lj.cpp b/descriptors/d2/d2_lj.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b0523480ce34325794b55aea3e7983028c2571c5
--- /dev/null
+++ b/descriptors/d2/d2_lj.cpp
@@ -0,0 +1,61 @@
+#include "d2_lj.h"
+
+Registry<D2_Base,Config&>::Register<D2_LJ> D2_LJ_1( "D2_LJ" );
+
+D2_LJ::D2_LJ(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{
+ if (!c.get<bool>("INIT2B")) {
+ s=0;
+ }
+}
+
+void D2_LJ::calc_aed(
+ const double,
+ const double rij_sq,
+ const double fc_ij,
+ aed_rtype aed)
+{
+ double r2_inv = 1.0/rij_sq;
+ double r6_inv = r2_inv*r2_inv*r2_inv;
+
+ aed(fidx) -= r6_inv*fc_ij;
+ aed(fidx+1) += r6_inv*r6_inv*fc_ij;
+}
+void D2_LJ::calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij)
+{
+
+ double r2_inv = 1.0/rij_sq;
+ double r6_inv = r2_inv*r2_inv*r2_inv;
+
+ fd_ij(fidx,0) = 6.0*r6_inv*r2_inv*rij*fc_ij - fcp_ij*r6_inv;
+ fd_ij(fidx+1,0) = -12.0*r6_inv*r6_inv*r2_inv*rij*fc_ij + fcp_ij*r6_inv*r6_inv;
+}
+void D2_LJ::calc_all(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij)
+{
+ double r2_inv = 1.0/rij_sq;
+ double r6_inv = r2_inv*r2_inv*r2_inv;
+
+ aed(fidx) -= r6_inv*fc_ij;
+ aed(fidx+1) += r6_inv*r6_inv*fc_ij;
+
+ fd_ij(fidx,0) = 6.0*r6_inv*r2_inv*rij*fc_ij - fcp_ij*r6_inv;
+ fd_ij(fidx+1,0) = -12.0*r6_inv*r6_inv*r2_inv*rij*fc_ij + fcp_ij*r6_inv*r6_inv;
+}
+size_t D2_LJ::size() {
+ return s;
+}
+std::string D2_LJ::label() {
+ return lab;
+}
diff --git a/descriptors/d2/d2_lj.h b/descriptors/d2/d2_lj.h
new file mode 100644
index 0000000000000000000000000000000000000000..a1a3f8dd6d3dfed2d4eb41924e30b7248dd5875c
--- /dev/null
+++ b/descriptors/d2/d2_lj.h
@@ -0,0 +1,59 @@
+#ifndef D2_LJ_H
+#define D2_LJ_H
+#include "d2_base.h"
+/** \anchor D2_LJ \brief Standard Lennard - Jones descriptor
+ *
+ * \f[
+ * V_i = \sum_{j \neq i} 4 \epsilon \Bigg(\Big(\frac{\sigma}{r_{ij}}\Big)^{12} - \Big(\frac{\sigma}{r_{ij}}\Big)^6\Bigg) f_c(r_{ij})
+ * \f]
+ *
+ * or equivalently:
+ *
+ * \f[
+ * V_i = \sum_{j \neq i} \frac{C_{12}}{r_{ij}^{12}} - \frac{C_6}{r_{ij}^6} f_c(r_{ij})
+ * \f]
+ *
+ * Note that machined learned coefficients \f$C_6\f$ and \f$C_{12}\f$
+ * corresponds to \f$\sigma\f$ and \f$\epsilon\f$ through the following relation:
+ *
+ * \f[ \sigma = \Big(\frac{C_{12}}{C_6}\Big)^{1/6} \f]
+ * \f[ \epsilon = \frac{1}{4} \frac{C_6^2}{C_{12}} w(Z) \f]
+ * where \f$w(Z)\f$ is a species depended weight factor (default is an atomic number).
+ *
+ * The machine learned \f$\sigma\f$ and \f$\epsilon\f$ only make sense
+ * (say to compare with the literature ones) when \ref BIAS false
+ * and \ref NORM false and system in monatomic.
+ * It is ok thought to set them to true it's just that
+ * numerical values will be different.
+ *
+ * Required Config Key: \ref INIT2B
+ */
+class D2_LJ : public D2_Base {
+ private:
+ size_t s=2;
+ std::string lab="D2_LJ";
+ int verbose;
+ public:
+ D2_LJ(Config &config);
+ void calc_aed(
+ const double,
+ const double rij_sq,
+ const double fc_ij,
+ aed_rtype aed);
+ void calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij);
+ void calc_all(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d2/d2_mie.cpp b/descriptors/d2/d2_mie.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d17366417d8a3e677120cd330588822189e8cc2f
--- /dev/null
+++ b/descriptors/d2/d2_mie.cpp
@@ -0,0 +1,72 @@
+#include "d2_mie.h"
+
+Registry<D2_Base,Config&>::Register<D2_MIE> D2_MIE_1( "D2_MIE" );
+
+D2_MIE::D2_MIE(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{
+ if (!c.get<bool>("INIT2B")) {
+ s=0;
+ return;
+ }
+ get_grid(c,"SGRID2B",etas);
+ if (etas.size()!=2) {
+ throw std::runtime_error("Number of elements in SGRID2B is != 2\n\
+ Mie descriptor requires SGRID2B with two positive elements.\n");
+ }
+ n=etas[1];
+ m=etas[0];
+ if (n<0 || m<0) {
+ throw std::runtime_error("Both Mie exponents must by positive\n");
+ }
+}
+
+void D2_MIE::calc_aed(
+ const double rij,
+ const double,
+ const double fc_ij,
+ aed_rtype aed)
+{
+ double rn_inv = 1.0/pow(rij,n);
+ double rm_inv = 1.0/pow(rij,m);
+
+ aed(fidx) -= rn_inv*fc_ij;
+ aed(fidx+1) += rm_inv*fc_ij;
+}
+void D2_MIE::calc_dXijdri(
+ const double rij,
+ const double,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij)
+{
+
+ fd_ij(fidx,0) = n/pow(rij,n+1)*fc_ij - fcp_ij/pow(rij,n);
+ fd_ij(fidx+1,0) = -m/pow(rij,m+1)*fc_ij + fcp_ij/pow(rij,m);
+
+}
+void D2_MIE::calc_all(
+ const double rij,
+ const double,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij)
+{
+
+ double rn_inv = 1.0/pow(rij,n);
+ double rm_inv = 1.0/pow(rij,m);
+
+ aed(fidx) -= rn_inv*fc_ij;
+ aed(fidx+1) += rm_inv*fc_ij;
+
+ fd_ij(fidx,0) = n*rn_inv/rij*fc_ij - fcp_ij*rn_inv;
+ fd_ij(fidx+1,0) = -m*rm_inv/rij*fc_ij + fcp_ij*rm_inv;
+
+}
+size_t D2_MIE::size() {
+ return s;
+}
+std::string D2_MIE::label() {
+ return lab;
+}
diff --git a/descriptors/d2/d2_mie.h b/descriptors/d2/d2_mie.h
new file mode 100644
index 0000000000000000000000000000000000000000..054fb221acbefb0774bac9a41e163eae4f27f160
--- /dev/null
+++ b/descriptors/d2/d2_mie.h
@@ -0,0 +1,54 @@
+#ifndef D2_MIE_H
+#define D2_MIE_H
+#include "d2_base.h"
+/** \brief Mie descriptor
+ *
+ * \f[
+ * V_i = \sum_{j \neq i} C \epsilon \Bigg(\Big(\frac{\sigma}{r_{ij}}\Big)^{n} - \Big(\frac{\sigma}{r_{ij}}\Big)^m\Bigg)
+ * \f]
+ *
+ * where
+ * \f[
+ * C=\frac{n}{n-m}\Big( \frac{n}{m} \Big)^{\frac{m}{n-m}}
+ * \f]
+ *
+ *
+ * Any cutoff can be used
+ *
+ * Required Config Key: \ref INIT2B \ref SGRID2B
+ *
+ * e.g. SGID2B 12 6
+ *
+ * will result in Lennard-Jones type descriptor
+ */
+class D2_MIE : public D2_Base {
+ private:
+ size_t s=2;
+ std::string lab="D2_MIE";
+ double m,n;
+ v_type etas;
+ int verbose;
+ public:
+ D2_MIE(Config &config);
+ void calc_aed(
+ const double,
+ const double rij_sq,
+ const double fc_ij,
+ aed_rtype aed);
+ void calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ fd_type &fd_ij);
+ void calc_all(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const double fcp_ij,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d3/d3_all.h b/descriptors/d3/d3_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..881cd74cb2651b737fa33d0e0a945b4ca34f84bf
--- /dev/null
+++ b/descriptors/d3/d3_all.h
@@ -0,0 +1,2 @@
+#include "d3_base.h"
+#include "d3_dummy.h"
diff --git a/descriptors/d3/d3_base.cpp b/descriptors/d3/d3_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f8c72a35e8a7989879412ea43aed0cafa944c623
--- /dev/null
+++ b/descriptors/d3/d3_base.cpp
@@ -0,0 +1,3 @@
+#include "d3_base.h"
+//template struct Registry<D3_Base, Config &>;
+
diff --git a/descriptors/d3/d3_base.h b/descriptors/d3/d3_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..281932a89db29a8ce063d0c4929ff66f04c6f0cf
--- /dev/null
+++ b/descriptors/d3/d3_base.h
@@ -0,0 +1,40 @@
+#ifndef D3_BASE_H
+#define D3_BASE_H
+
+#include "../d_base.h"
+
+class D3_Base: public D_Base {
+ public:
+ virtual ~D3_Base() {};
+
+ virtual void calc_aed(
+ const size_t fidx,
+ const double rij,
+ const double rik,
+ const double fc_ij,
+ const double fc_ik,
+ aed_rtype aed)=0;
+ virtual void calc_fd(
+ const size_t fidx,
+ const double rij,
+ const double rik,
+ const double fc_ij,
+ const double fc_ik,
+ const double fcp_ij,
+ const double fcp_ik,
+ fd_type &fd_ij)=0;
+ virtual void calc_all(
+ const size_t fidx,
+ const double rij,
+ const double rik,
+ const double fc_ij,
+ const double fc_ik,
+ const double fcp_ij,
+ const double fcp_ik,
+ aed_rtype aed,
+ fd_type &fd_ij)=0;
+ virtual size_t size()=0;
+ virtual std::string label()=0;
+};
+template<> inline Registry<D3_Base,Config&>::Map Registry<D3_Base,Config&>::registry{};
+#endif
diff --git a/descriptors/d3/d3_dummy.cpp b/descriptors/d3/d3_dummy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c145bfd552d087e758dbe6ae6d6fe11d5982b66f
--- /dev/null
+++ b/descriptors/d3/d3_dummy.cpp
@@ -0,0 +1,35 @@
+#include "d3_dummy.h"
+
+Registry<D3_Base,Config&>::Register<D3_Dummy> D3_Dummy_1( "D3_Dummy" );
+
+D3_Dummy::D3_Dummy() {}
+D3_Dummy::D3_Dummy(Config &) {}
+
+void D3_Dummy::calc_aed(
+ const size_t,
+ const double,
+ const double,
+ const double,
+ const double,
+ aed_rtype ) {}
+void D3_Dummy::calc_fd(
+ const size_t,
+ const double,
+ const double,
+ const double,
+ const double,
+ const double,
+ const double,
+ fd_type &) {}
+void D3_Dummy::calc_all(
+ const size_t,
+ const double,
+ const double,
+ const double,
+ const double,
+ const double,
+ const double,
+ aed_rtype ,
+ fd_type &) {}
+size_t D3_Dummy::size() { return s;}
+std::string D3_Dummy::label() {return lab;}
diff --git a/descriptors/d3/d3_dummy.h b/descriptors/d3/d3_dummy.h
new file mode 100644
index 0000000000000000000000000000000000000000..616b08a021e7e6159589baf8c9d1355b1f66cb37
--- /dev/null
+++ b/descriptors/d3/d3_dummy.h
@@ -0,0 +1,40 @@
+#ifndef D3_DUMMY_H
+#define D3_DUMMY_H
+#include "d3_base.h"
+class D3_Dummy: public D3_Base {
+ private:
+ size_t s=0;
+ std::string lab="D3_Dummy";
+ public:
+ D3_Dummy();
+ D3_Dummy(Config&);
+ void calc_aed(
+ const size_t fidx,
+ const double rij,
+ const double rik,
+ const double fc_ij,
+ const double fc_ik,
+ aed_rtype aed);
+ void calc_fd(
+ const size_t fidx,
+ const double rij,
+ const double rik,
+ const double fc_ij,
+ const double fc_ik,
+ const double fcp_ij,
+ const double fcp_ik,
+ fd_type &fd_ij);
+ void calc_all(
+ const size_t fidx,
+ const double rij,
+ const double rik,
+ const double fc_ij,
+ const double fc_ik,
+ const double fcp_ij,
+ const double fcp_ik,
+ aed_rtype aed,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+};
+#endif
diff --git a/descriptors/d_all.h b/descriptors/d_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..f877cc369427025c2cc11561b252217058885296
--- /dev/null
+++ b/descriptors/d_all.h
@@ -0,0 +1,3 @@
+#include "d2/d2_all.h"
+#include "d3/d3_all.h"
+#include "dm/dm_all.h"
diff --git a/descriptors/d_base.h b/descriptors/d_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ac2101226451074a6a553ddbfc3edde6c03015e
--- /dev/null
+++ b/descriptors/d_base.h
@@ -0,0 +1,57 @@
+#ifndef D_BASE_H
+#define D_BASE_H
+
+#include "../../CORE/config/config.h"
+#include "../../CORE/typedefs.h"
+#include "../../LIBS/Eigen/Dense"
+#include "../../CORE/registry.h"
+
+/** \brief Base class for all descriptor types.
+ *
+ */
+class D_Base {
+ public:
+ static void get_grid(const Config &c, std::string key, v_type &v) {
+ if (c.get<double>(key) < 0) {
+ // generate grid automatically if first element
+ // is<zero and int. If it is double than assume that
+ // grid is provided manually and min value is just smaller
+ // than zero.
+ double int_part;
+ if (std::modf ( c.get<double>(key), &int_part ) == 0.0) {
+ // automatic generation
+ size_t cg = abs(c.get<int>(key,1));
+ double cgmin = c.get<double>(key,2);
+ double cgmax = c.get<double>(key,3);
+ if (c.get<int>(key) == -1) {
+ v = linspace(cgmin, cgmax, cg);
+ }
+ else if (c.get<int>(key) == -2) {
+ v = logspace(cgmin, cgmax, cg, exp(1));
+ }
+ else {
+ throw std::runtime_error(key+" algorithm not supported\n");
+ }
+ }
+ else {
+ v.resize(c.size(key));
+ c.get<v_type>(key,v);
+ }
+ }
+ else {
+ v.resize(c.size(key));
+ c.get<v_type>(key,v);
+ }
+ }
+ int verbose;
+ virtual ~D_Base() {};
+
+ /** \brief Return dimension of the descriptor.
+ */
+ virtual size_t size()=0;
+
+ /** \brief Return label of this descriptor.
+ */
+ virtual std::string label()=0;
+};
+#endif
diff --git a/descriptors/d_basis_functions.h b/descriptors/d_basis_functions.h
new file mode 100644
index 0000000000000000000000000000000000000000..c464210b487869b3f552417839ce79f5857f35fc
--- /dev/null
+++ b/descriptors/d_basis_functions.h
@@ -0,0 +1,101 @@
+#ifndef D_BASIS_F_H
+#define D_BASIS_F_H
+#include <math.h>
+
+// GAUSSIANS
+inline double G(double r, double eta, double miu) {
+ return exp(-eta*(r-miu)*(r-miu));
+}
+inline double G(double r, double eta, double miu, double f) {
+ return exp(-eta*(r-miu)*(r-miu))*f;
+}
+inline double dG(double r, double eta, double miu) {
+ return -2.0*eta*(r-miu)*exp(-eta*(r-miu)*(r-miu));
+}
+inline double dG(double r, double eta, double miu, double f, double fp) {
+ return exp(-eta*(r-miu)*(r-miu) )*(-2.0*f*eta*(r-miu)+fp);
+}
+
+// BLIPS
+inline double B(double x_c)
+ // x_c = eta*(rij-r_s)
+ // return B(x_c)
+{
+ if (-2.0 < x_c && x_c <= -1.0) {
+ double b = 2.0 + x_c;
+ return 0.25 * b*b*b;
+ }
+ else if (-1.0 < x_c && x_c <= 1.0) {
+ double x2 = x_c*x_c;
+ double x3 = x2*x_c;
+ return 1.0 - 1.5*x2 + 0.75*std::fabs(x3);
+ }
+ else if (1.0 < x_c && x_c < 2.0) {
+ double b = 2.0 - x_c;
+ return 0.25 * b*b*b;
+ }
+ else {
+ return 0.0;
+ }
+}
+inline double B(double x_c, double f)
+ // x_c = eta*(rij-r_s)
+ // return B(x_c)*f
+{
+ if (-2.0 < x_c && x_c <= -1.0) {
+ double b = 2.0 + x_c;
+ return 0.25 * b*b*b *f;
+ }
+ else if (-1.0 < x_c && x_c <= 1.0) {
+ double x2 = x_c*x_c;
+ double x3 = x2*x_c;
+ return f*(1.0 - 1.5*x2 + 0.75*std::fabs(x3));
+ }
+ else if (1.0 < x_c && x_c < 2.0) {
+ double b = 2.0 - x_c;
+ return f*(0.25 * b*b*b);
+ }
+ else {
+ return 0.0;
+ }
+}
+inline double dB(double x_c, double eta)
+ // def: x_c = eta*(rij-r_s)
+ // return d/dr_ij B(x_c) = eta*d/dx_c B(x_c)
+{
+ if (-2.0 < x_c && x_c <= -1.0) {
+ double d = 2.0+x_c;
+ return 0.75*eta*d*d;
+ }
+ else if (-1.0 < x_c && x_c <= 1.0) {
+ return -3.0*eta*x_c + 2.25*eta*x_c*std::fabs(x_c);
+ }
+ else if (1.0 < x_c && x_c < 2.0) {
+ double d = 2.0-x_c;
+ return -0.75*eta*d*d;
+ }
+ else {
+ return 0.0;
+ }
+}
+inline double dB(double x_c, double eta, double f, double fp)
+ // def: x_c = eta*(rij-r_s)
+ // return d/dr_ij f(r_ij)*B(x_c)
+{
+ if (-2.0 < x_c && x_c <= -1.0) {
+ double d = 2.0+x_c;
+ return f*(0.75*eta*d*d)+fp*(0.25 * d*d*d);
+ }
+ else if (-1.0 < x_c && x_c <= 1.0) {
+ return f*(-3.0*eta*x_c + 2.25*eta*x_c*std::fabs(x_c))
+ +fp*(1.0 - 1.5*x_c*x_c + 0.75*std::fabs(x_c*x_c*x_c));
+ }
+ else if (1.0 < x_c && x_c < 2.0) {
+ double d = 2.0-x_c;
+ return f*(-0.75*eta*d*d)+fp*(0.25 * d*d*d);
+ }
+ else {
+ return 0.0;
+ }
+}
+#endif
diff --git a/descriptors/dm/dm_all.h b/descriptors/dm/dm_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..5f14169c7a5fedff1ecef0869be62f09a6eacd8a
--- /dev/null
+++ b/descriptors/dm/dm_all.h
@@ -0,0 +1,6 @@
+#include "dm_base.h"
+#include "dm_dummy.h"
+#include "dm_eam.h"
+#include "dm_ead.h"
+#include "dm_blip.h"
+#include "dm_mEAD.h"
diff --git a/descriptors/dm/dm_base.cpp b/descriptors/dm/dm_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0099681b1582d07b610298e9f0cd79bfa451029b
--- /dev/null
+++ b/descriptors/dm/dm_base.cpp
@@ -0,0 +1,2 @@
+#include "dm_base.h"
+
diff --git a/descriptors/dm/dm_base.h b/descriptors/dm/dm_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..1912e4be8d4c463cbb018073211dad5c474b9eaf
--- /dev/null
+++ b/descriptors/dm/dm_base.h
@@ -0,0 +1,90 @@
+#ifndef DM_BASE_H
+#define DM_BASE_H
+
+#include <string>
+#include "../d_base.h"
+
+/** \brief Base class for all many-body type descriptors.
+ *
+ * All many-body descriptors **must** inherit this class.
+ */
+class DM_Base: public D_Base {
+ public:
+ size_t fidx=0; // first index in aed and fd arrays for this descriptor
+ virtual ~DM_Base() {};
+
+ /** \brief Calculate \ref AED
+ *
+ * Calculate Atomic Energy Descriptor for the atom local environment.
+ */
+ virtual void calc_aed(
+ rho_rtype rho,
+ aed_rtype aed)=0;
+
+ /** \brief Calculate \ref FD
+ *
+ * Calculate Force Descriptor between to atoms.
+ *
+ * This method works for half NN lists
+ * and linear models.
+ *
+ * \f[
+ * \mathbf{fd}_{ij} \mathrel{+}=
+ * \frac{\partial \mathbf{X}_{ij}}{\partial \mathbf{r}_i} +
+ * \frac{\partial \mathbf{X}_{ji}}{\partial \mathbf{r}_i}
+ * \f]
+ */
+ virtual int calc_dXijdri_dXjidri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj)=0;
+
+ /** \brief Calculate \ref FD
+ *
+ * Calculate Force Descriptor between to atoms.
+ *
+ * This method works for full NN lists
+ * and all models.
+ *
+ * \f[
+ * \mathbf{fd}_{ij} \mathrel{+}=
+ * \frac{\partial \mathbf{X}_{ij}}{\partial \mathbf{r}_i}
+ * \f]
+ */
+ virtual int calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij)=0;
+
+ /** \brief Resize arrays for a density and F' */
+ virtual void init_rhoi(rho_rtype rhoi)=0;
+
+ /** \brief Calculate density */
+ virtual void calc_rho(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rho)=0;
+
+ /** \brief Return size of the density array */
+ virtual size_t rhoi_size()=0;
+
+ /** \brief Return size of the derivative of the embedding energy array */
+ virtual size_t rhoip_size()=0;
+ virtual size_t size()=0;
+ virtual std::string label()=0;
+};
+template<> inline Registry<DM_Base,Config&>::Map Registry<DM_Base,Config&>::registry{};
+#endif
diff --git a/descriptors/dm/dm_blip.cpp b/descriptors/dm/dm_blip.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0a23ea92e4c3c98e5a313643091f288c56d2bb2a
--- /dev/null
+++ b/descriptors/dm/dm_blip.cpp
@@ -0,0 +1,299 @@
+#include "dm_blip.h"
+#include "../d_basis_functions.h"
+#include <vector>
+#include <iomanip>
+
+Registry<DM_Base,Config&>::Register<DM_Blip> DM_Blip_1( "DM_Blip" );
+
+DM_Blip::DM_Blip(Config &c):
+ verbose(c.get<int>("VERBOSE")),
+ rc(c.get<double>("RCUTMB"))
+{
+
+ if (!c.get<bool>("INITMB")) return;
+
+ get_grid(c,"CGRIDMB",cgrid);
+ get_grid(c,"SGRIDMB",sgrid);
+ if (cgrid.size()!=sgrid.size()) {
+
+ throw std::runtime_error("SGRID2B and CGRID2B arrays differ in size.\n");
+ }
+
+ Lmax = c.get<int>("AGRIDMB");
+ s=sgrid.size()*(Lmax+1);
+ rhoisize = gen_atomic_orbitals(Lmax);
+ rhoisize *= sgrid.size();
+
+ if (verbose) {
+ std::cout << std::endl;
+ std::cout << "SGRID (SGRIDMB): " << sgrid.size() << std::endl;
+ for (auto e:sgrid) std::cout << e << " ";
+ std::cout << std::endl;
+
+ std::cout << "CGRID (CGRIDMB): " << cgrid.size() << std::endl;
+ for (auto L:cgrid) std::cout << L << " ";
+ std::cout << std::endl;
+ std::cout << "rhoisize: " << rhoisize << std::endl;
+ }
+}
+
+void DM_Blip::calc_aed(
+ rho_rtype rhoi,
+ aed_rtype aed)
+{
+ size_t Lshift = 0;
+ size_t ii=0;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const double f = fac[L][o];
+ for (size_t c=0; c<cgrid.size(); c++) {
+ aed(c+Lshift+fidx) += f*rhoi(ii)*rhoi(ii);
+ rhoi(ii+rhoisize) = 2.0*f*rhoi(ii);
+ ii++;
+ }
+ }
+ Lshift += sgrid.size();
+ }
+}
+int DM_Blip::calc_dXijdri_dXjidri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj)
+{
+ const double x = vec_ij[0];
+ const double y = vec_ij[1];
+ const double z = vec_ij[2];
+
+ size_t Lshift = 0;
+
+ size_t ii=rhoisize;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const int lx = orbitals[L][o][0];
+ const int ly = orbitals[L][o][1];
+ const int lz = orbitals[L][o][2];
+
+ const double powxlxij = my_pow(x,lx);
+ const double powylyij = my_pow(y,ly);
+ const double powzlzij = my_pow(z,lz);
+
+ const double powxlxji = my_pow(-x,lx);
+ const double powylyji = my_pow(-y,ly);
+ const double powzlzji = my_pow(-z,lz);
+
+ double txij=0.0,tyij=0.0,tzij=0.0;
+ double txji=0.0,tyji=0.0,tzji=0.0;
+ if (lx!=0) {
+ txij = lx*my_pow(x,lx-1)*powylyij*powzlzij;
+ txji = lx*my_pow(-x,lx-1)*powylyji*powzlzji;
+ }
+ if (ly!=0) {
+ tyij = ly*my_pow(y,ly-1)*powxlxij*powzlzij;
+ tyji = ly*my_pow(-y,ly-1)*powxlxji*powzlzji;
+ }
+ if (lz!=0) {
+ tzij = lz*my_pow(z,lz-1)*powylyij*powxlxij;
+ tzji = lz*my_pow(-z,lz-1)*powylyji*powxlxji;
+ }
+ const double pqrij = powxlxij*powylyij*powzlzij;
+ const double pqrji = powxlxji*powylyji*powzlzji;
+
+ for (size_t c=0; c<cgrid.size(); c++) {
+ double blip = B(sgrid[c]*(rij-cgrid[c]),fc_ij);
+ double dblip = dB(sgrid[c]*(rij-cgrid[c]),sgrid[c],fc_ij,fcp_ij);
+ const double term2ijx=pqrij*dblip*x/rij;
+ const double term2ijy=pqrij*dblip*y/rij;
+ const double term2ijz=pqrij*dblip*z/rij;
+
+ const double term2jix=-pqrji*dblip*x/rij;
+ const double term2jiy=-pqrji*dblip*y/rij;
+ const double term2jiz=-pqrji*dblip*z/rij;
+
+ const double term1ijx = blip*txij;
+ const double term1ijy = blip*tyij;
+ const double term1ijz = blip*tzij;
+
+ const double term1jix = blip*txji;
+ const double term1jiy = blip*tyji;
+ const double term1jiz = blip*tzji;
+
+ fd_ij(c+Lshift+fidx,0) +=
+ rhoi(ii)*(term1ijx + term2ijx)*wj
+ -rhoj(ii)*(term1jix + term2jix)*wi
+ ;
+
+ fd_ij(c+Lshift+fidx,1) +=
+ rhoi(ii)*(term1ijy + term2ijy)*wj
+ -rhoj(ii)*(term1jiy + term2jiy)*wi
+ ;
+
+ fd_ij(c+Lshift+fidx,2) +=
+ rhoi(ii)*(term1ijz + term2ijz)*wj
+ -rhoj(ii)*(term1jiz + term2jiz)*wi
+ ;
+ ii++;
+ }
+ }
+ Lshift+=sgrid.size();
+ }
+
+ return 1;
+
+}
+int DM_Blip::calc_dXijdri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij)
+{
+ const double x = vec_ij[0];
+ const double y = vec_ij[1];
+ const double z = vec_ij[2];
+
+ size_t Lshift = 0;
+
+ size_t ii=rhoisize;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const int lx = orbitals[L][o][0];
+ const int ly = orbitals[L][o][1];
+ const int lz = orbitals[L][o][2];
+
+ const double powxlxij = my_pow(x,lx);
+ const double powylyij = my_pow(y,ly);
+ const double powzlzij = my_pow(z,lz);
+
+ double txij=0.0,tyij=0.0,tzij=0.0;
+ if (lx!=0) {
+ txij = lx*my_pow(x,lx-1)*powylyij*powzlzij;
+ }
+ if (ly!=0) {
+ tyij = ly*my_pow(y,ly-1)*powxlxij*powzlzij;
+ }
+ if (lz!=0) {
+ tzij = lz*my_pow(z,lz-1)*powylyij*powxlxij;
+ }
+ const double pqrij = powxlxij*powylyij*powzlzij;
+
+ for (size_t c=0; c<cgrid.size(); c++) {
+ double blip = B(sgrid[c]*(rij-cgrid[c]),fc_ij);
+ double dblip = dB(sgrid[c]*(rij-cgrid[c]),sgrid[c],fc_ij,fcp_ij);
+
+ const double term2ijx=pqrij*dblip*x/rij;
+ const double term2ijy=pqrij*dblip*y/rij;
+ const double term2ijz=pqrij*dblip*z/rij;
+
+ const double term1ijx = blip*txij;
+ const double term1ijy = blip*tyij;
+ const double term1ijz = blip*tzij;
+
+ fd_ij(c+Lshift+fidx,0) +=
+ rhoi(ii)*(term1ijx + term2ijx)
+ ;
+
+ fd_ij(c+Lshift+fidx,1) +=
+ rhoi(ii)*(term1ijy + term2ijy)
+ ;
+
+ fd_ij(c+Lshift+fidx,2) +=
+ rhoi(ii)*(term1ijz + term2ijz)
+ ;
+ ii++;
+ }
+ }
+ Lshift+=sgrid.size();
+ }
+
+ return 1;
+
+}
+size_t DM_Blip::size() {
+ return s;
+}
+std::string DM_Blip::label() {
+ return lab;
+}
+
+void DM_Blip::init_rhoi(rho_rtype rhoi)
+{
+ rhoi.resize(2*rhoisize);
+ rhoi.setZero();
+}
+void DM_Blip::calc_rho(
+ const double rij,
+ const double,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rhoi)
+{
+ size_t ii=0;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const size_t lx = orbitals[L][o][0];
+ const size_t ly = orbitals[L][o][1];
+ const size_t lz = orbitals[L][o][2];
+ double t = (my_pow(vec_ij[0],lx)*my_pow(vec_ij[1],ly)*my_pow(vec_ij[2],lz));
+ for (size_t c=0; c<cgrid.size(); c++) {
+ rhoi(ii++) += B(sgrid[c]*(rij-cgrid[c]),fc_ij)*t;
+ }
+ }
+ }
+}
+
+size_t DM_Blip::gen_atomic_orbitals(int Lmax) {
+ orbitals.resize(Lmax+1);
+ fac.resize(Lmax+1);
+
+ size_t count = 0;
+ for (int L = 0; L <= Lmax; ++L) {
+ for (int lx = 0; lx <= L; ++lx) {
+ for (int ly = 0; ly <= L; ++ly) {
+ for (int lz = 0; lz <= L; ++lz) {
+ if (lx+ly+lz == L) {
+ std::vector<int> o = {lx, ly, lz};
+ orbitals[L].push_back(o);
+ const double f = fact(L)/(fact(lx)*fact(ly)*fact(lz))/my_pow(4*rc,L);
+ fac[L].push_back(f);
+ count++;
+ }
+ }
+ }
+ }
+ }
+ return count;
+}
+double DM_Blip::fact(int n)
+{
+ double f=1.0;
+ while (n) {
+ f *= n;
+ --n;
+ }
+ return f;
+}
+double DM_Blip::my_pow(double x, int n){
+ double r = 1.0;
+
+ while(n > 0){
+ r *= x;
+ --n;
+ }
+
+ return r;
+}
+size_t DM_Blip::rhoi_size() {
+ return rhoisize;
+}
+size_t DM_Blip::rhoip_size() {
+ return rhoisize;
+}
diff --git a/descriptors/dm/dm_blip.h b/descriptors/dm/dm_blip.h
new file mode 100644
index 0000000000000000000000000000000000000000..f11e443cb746f0954bfa37c87de8e00c8e5120cd
--- /dev/null
+++ b/descriptors/dm/dm_blip.h
@@ -0,0 +1,90 @@
+#ifndef DM_BLIP_H
+#define DM_BLIP_H
+#include "dm_base.h"
+
+/** Blip Many Body Descriptor.
+ *
+ * \f[
+ * V_i^{L,\eta,r_s} = \sum_{l_x,l_y,l_z}^{l_x+l_y+l_z=L} \frac{L!}{l_x!l_y!l_z!}
+ * \Big( \rho_i^{\eta,r_s,l_x,l_y,l_z} \Big)^2
+ * \f]
+ *
+ * where density \f$ \rho \f$ is calculated using modified Gaussian Type Orbitals
+ * (expansion in the Blip basis instead of usual Gaussians):
+ * \f[
+ * \rho_i^{\eta,r_s,l_x,l_y,l_z} = \sum_{j \neq i} x_{ij}^{l_x}y_{ij}^{l_y}z_{ij}^{l_z}
+ * \B{\Big(-\eta(r_{ij}-r_s)^2\Big)}f_c(r_{ij})
+ * \f]
+ *
+ * \ref CGRIDMB parameters control position \f$ r_s \f$ of the gaussian basis function.
+ *
+ * \ref SGRIDMB parameters control width \f$ \eta \f$ of the gaussian basis function.
+ *
+ * \ref AGRIDMB parameter specify maximum value for the angular momentum \f$L_{max}\f$.
+ *
+ * e.g. \f$L_{max}=2\f$ will calculate descriptors with \f$ L=0,1,2 \f$ (s,p,d orbitals).
+ *
+ * More information about this descriptor:
+ *
+ * <div class="csl-entry">Zhang, Y., Hu, C.,Jiang, B.
+ * (2019). Embedded atom neural network potentials: efficient and accurate
+ * machine learning with a physically inspired representation.
+ * <i>Journal of Physical Chemistry Letters</i>, <i>10</i>(17),
+ * 4962–4967. https://doi.org/10.1021/acs.jpclett.9b02037</div>
+ *
+ * Required Config keys:
+ * \ref INITMB \ref CGRIDMB \ref SGRIDMB \ref AGRIDMB
+ */
+class DM_Blip: public DM_Base {
+ private:
+ int Lmax=0;
+ size_t rhoisize=0;
+ size_t s=0;
+ std::string lab="DM_Blip";
+ v_type sgrid;
+ v_type cgrid;
+ std::vector<std::vector<std::vector<int>>> orbitals;
+ std::vector<std::vector<double>> fac;
+ size_t gen_atomic_orbitals(int L);
+ double fact(int n);
+ double my_pow(double,int);
+ int verbose;
+ double rc;
+
+ public:
+ DM_Blip(Config&);
+ void calc_aed(
+ rho_rtype rho,
+ aed_rtype aed);
+ int calc_dXijdri_dXjidri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj);
+ int calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+ void init_rhoi(rho_rtype rhoi);
+ void calc_rho(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rho);
+ size_t rhoi_size();
+ size_t rhoip_size();
+};
+#endif
diff --git a/descriptors/dm/dm_dummy.cpp b/descriptors/dm/dm_dummy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eff22e7581d0eee12deda7a35a385a67ed3e3f27
--- /dev/null
+++ b/descriptors/dm/dm_dummy.cpp
@@ -0,0 +1,40 @@
+#include "dm_dummy.h"
+
+Registry<DM_Base,Config&>::Register<DM_Dummy> DM_Dummy_1( "DM_Dummy" );
+
+DM_Dummy::DM_Dummy() {}
+DM_Dummy::DM_Dummy(Config&) {}
+void DM_Dummy::calc_aed(
+ rho_rtype,
+ aed_rtype ) {}
+int DM_Dummy::calc_dXijdri_dXjidri(
+ const double,
+ const double,
+ const Eigen::Vector3d &,
+ const double,
+ const double,
+ rho_rtype,
+ rho_rtype,
+ fd_type &,
+ const double,
+ const double) { return 0; }
+int DM_Dummy::calc_dXijdri(
+ const double,
+ const double,
+ const Eigen::Vector3d &,
+ const double,
+ const double,
+ rho_rtype,
+ fd_type &) { return 0; }
+size_t DM_Dummy::size() { return s;}
+std::string DM_Dummy::label() {return lab;}
+void DM_Dummy::init_rhoi(rho_rtype) {}
+void DM_Dummy::calc_rho(
+ const double,
+ const double,
+ const double,
+ const Eigen::Vector3d &,
+ rho_rtype) {}
+size_t DM_Dummy::rhoi_size() { return 0; }
+size_t DM_Dummy::rhoip_size() { return 0; }
+
diff --git a/descriptors/dm/dm_dummy.h b/descriptors/dm/dm_dummy.h
new file mode 100644
index 0000000000000000000000000000000000000000..9996862bef89b185e7bc6cff5ff29e71de5e4429
--- /dev/null
+++ b/descriptors/dm/dm_dummy.h
@@ -0,0 +1,53 @@
+#ifndef DM_DUMMY_H
+#define DM_DUMMY_H
+#include "dm_base.h"
+
+/**
+ * Dummy many-body descriptor.
+ *
+ * Use it to satisfy DescriptorsCalc requirements in case
+ * when many-body descriptor is not required.
+ *
+ */
+class DM_Dummy: public DM_Base {
+ private:
+ size_t s=0;
+ std::string lab="DM_Dummy";
+ public:
+ DM_Dummy();
+ DM_Dummy(Config&);
+ void calc_aed(
+ rho_rtype rho,
+ aed_rtype aed);
+ int calc_dXijdri_dXjidri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj);
+ int calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+ void init_rhoi(rho_rtype rhoi);
+ void calc_rho(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rho);
+ size_t rhoi_size();
+ size_t rhoip_size();
+};
+#endif
diff --git a/descriptors/dm/dm_ead.cpp b/descriptors/dm/dm_ead.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..81045d9db7b7dbae8e7d29d68ca833ec92acfeb3
--- /dev/null
+++ b/descriptors/dm/dm_ead.cpp
@@ -0,0 +1,303 @@
+#include "dm_ead.h"
+#include "../d_basis_functions.h"
+#include <vector>
+#include <iomanip>
+
+Registry<DM_Base,Config&>::Register<DM_EAD> DM_EAD_1( "DM_EAD" );
+
+/* Index ii iterates rhoi array
+ * rhoi array is split 50/50 into density and derivative
+ * of the embedding function
+ */
+DM_EAD::DM_EAD(Config &config):
+ verbose(config.get<int>("VERBOSE")),
+ rc(config.get<double>("RCUTMB"))
+{
+
+ if (!config.get<bool>("INITMB")) return;
+
+ get_grid(config,"CGRIDMB",cgrid);
+ get_grid(config,"SGRIDMB",sgrid);
+
+ if (cgrid.size()!=sgrid.size()) {
+ throw std::runtime_error("SGRID2B and CGRID2B arrays differ in size.\n");
+ }
+
+ Lmax = config.get<int>("AGRIDMB");
+ s=sgrid.size()*(Lmax+1);
+ rhoisize = gen_atomic_orbitals(Lmax);
+ rhoisize *= sgrid.size();
+
+ if (verbose) {
+ std::cout << std::endl;
+ std::cout << "SGRID (SGRIDMB): " << sgrid.size() << std::endl;
+ for (auto e:sgrid) std::cout << e << " ";
+ std::cout << std::endl;
+
+ std::cout << "CGRID (CGRIDMB): " << cgrid.size() << std::endl;
+ for (auto L:cgrid) std::cout << L << " ";
+ std::cout << std::endl;
+ std::cout << "rhoisize: " << rhoisize << std::endl;
+ }
+}
+
+void DM_EAD::calc_aed(
+ rho_rtype rhoi,
+ aed_rtype aed)
+{
+ size_t Lshift = 0;
+ size_t ii=0;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const double f = fac[L][o];
+ for (size_t c=0; c<cgrid.size(); c++) {
+ aed(c+Lshift+fidx) += f*rhoi(ii)*rhoi(ii);
+ rhoi(ii+rhoisize) = 2.0*f*rhoi(ii);
+ ii++;
+ }
+ }
+ Lshift += sgrid.size();
+ }
+}
+int DM_EAD::calc_dXijdri_dXjidri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj)
+{
+ const double x = vec_ij[0];
+ const double y = vec_ij[1];
+ const double z = vec_ij[2];
+
+ size_t Lshift = 0;
+
+ size_t ii=rhoisize;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const int lx = orbitals[L][o][0];
+ const int ly = orbitals[L][o][1];
+ const int lz = orbitals[L][o][2];
+
+ const double powxlxij = my_pow(x,lx);
+ const double powylyij = my_pow(y,ly);
+ const double powzlzij = my_pow(z,lz);
+
+ const double powxlxji = my_pow(-x,lx);
+ const double powylyji = my_pow(-y,ly);
+ const double powzlzji = my_pow(-z,lz);
+
+ double txij=0.0,tyij=0.0,tzij=0.0;
+ double txji=0.0,tyji=0.0,tzji=0.0;
+ if (lx!=0) {
+ txij = lx*my_pow(x,lx-1)*powylyij*powzlzij;
+ txji = lx*my_pow(-x,lx-1)*powylyji*powzlzji;
+ }
+ if (ly!=0) {
+ tyij = ly*my_pow(y,ly-1)*powxlxij*powzlzij;
+ tyji = ly*my_pow(-y,ly-1)*powxlxji*powzlzji;
+ }
+ if (lz!=0) {
+ tzij = lz*my_pow(z,lz-1)*powylyij*powxlxij;
+ tzji = lz*my_pow(-z,lz-1)*powylyji*powxlxji;
+ }
+ const double pqrij = powxlxij*powylyij*powzlzij;
+ const double pqrji = powxlxji*powylyji*powzlzji;
+
+ for (size_t c=0; c<cgrid.size(); c++) {
+ double gauss = G(rij,sgrid[c],cgrid[c],fc_ij);
+ double dgauss = dG(rij,sgrid[c],cgrid[c],fc_ij,fcp_ij);
+
+ const double term2ijx=pqrij*dgauss*x/rij;
+ const double term2ijy=pqrij*dgauss*y/rij;
+ const double term2ijz=pqrij*dgauss*z/rij;
+
+ const double term2jix=-pqrji*dgauss*x/rij;
+ const double term2jiy=-pqrji*dgauss*y/rij;
+ const double term2jiz=-pqrji*dgauss*z/rij;
+
+ const double term1ijx = gauss*txij;
+ const double term1ijy = gauss*tyij;
+ const double term1ijz = gauss*tzij;
+
+ const double term1jix = gauss*txji;
+ const double term1jiy = gauss*tyji;
+ const double term1jiz = gauss*tzji;
+
+ fd_ij(c+Lshift+fidx,0) +=
+ rhoi(ii)*(term1ijx + term2ijx)*wj
+ -rhoj(ii)*(term1jix + term2jix)*wi
+ ;
+
+ fd_ij(c+Lshift+fidx,1) +=
+ rhoi(ii)*(term1ijy + term2ijy)*wj
+ -rhoj(ii)*(term1jiy + term2jiy)*wi
+ ;
+
+ fd_ij(c+Lshift+fidx,2) +=
+ rhoi(ii)*(term1ijz + term2ijz)*wj
+ -rhoj(ii)*(term1jiz + term2jiz)*wi
+ ;
+ ii++;
+ }
+ }
+ Lshift+=sgrid.size();
+ }
+
+ return 1;
+
+}
+int DM_EAD::calc_dXijdri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij)
+{
+ const double x = vec_ij[0];
+ const double y = vec_ij[1];
+ const double z = vec_ij[2];
+
+ size_t Lshift = 0;
+
+ size_t ii=rhoisize;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const int lx = orbitals[L][o][0];
+ const int ly = orbitals[L][o][1];
+ const int lz = orbitals[L][o][2];
+
+ const double powxlxij = my_pow(x,lx);
+ const double powylyij = my_pow(y,ly);
+ const double powzlzij = my_pow(z,lz);
+
+ double txij=0.0,tyij=0.0,tzij=0.0;
+ if (lx!=0) {
+ txij = lx*my_pow(x,lx-1)*powylyij*powzlzij;
+ }
+ if (ly!=0) {
+ tyij = ly*my_pow(y,ly-1)*powxlxij*powzlzij;
+ }
+ if (lz!=0) {
+ tzij = lz*my_pow(z,lz-1)*powylyij*powxlxij;
+ }
+ const double pqrij = powxlxij*powylyij*powzlzij;
+
+ for (size_t c=0; c<cgrid.size(); c++) {
+ double gauss = G(rij,sgrid[c],cgrid[c],fc_ij);
+ double dgauss = dG(rij,sgrid[c],cgrid[c],fc_ij,fcp_ij);
+
+ const double term2ijx=pqrij*dgauss*x/rij;
+ const double term2ijy=pqrij*dgauss*y/rij;
+ const double term2ijz=pqrij*dgauss*z/rij;
+
+ const double term1ijx = gauss*txij;
+ const double term1ijy = gauss*tyij;
+ const double term1ijz = gauss*tzij;
+
+ fd_ij(c+Lshift+fidx,0) +=
+ rhoi(ii)*(term1ijx + term2ijx)
+ ;
+
+ fd_ij(c+Lshift+fidx,1) +=
+ rhoi(ii)*(term1ijy + term2ijy)
+ ;
+
+ fd_ij(c+Lshift+fidx,2) +=
+ rhoi(ii)*(term1ijz + term2ijz)
+ ;
+ ii++;
+ }
+ }
+ Lshift+=sgrid.size();
+ }
+
+ return 1;
+
+}
+size_t DM_EAD::size() {
+ return s;
+}
+std::string DM_EAD::label() {
+ return lab;
+}
+
+void DM_EAD::init_rhoi(rho_rtype rhoi)
+{
+ rhoi.resize(2*rhoisize);
+ rhoi.setZero();
+}
+void DM_EAD::calc_rho(
+ const double rij,
+ const double,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rhoi)
+{
+ size_t ii=0;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const size_t lx = orbitals[L][o][0];
+ const size_t ly = orbitals[L][o][1];
+ const size_t lz = orbitals[L][o][2];
+ double t = (my_pow(vec_ij[0],lx)*my_pow(vec_ij[1],ly)*my_pow(vec_ij[2],lz));
+ for (size_t c=0; c<cgrid.size(); c++) {
+ rhoi(ii++) += G(rij,sgrid[c],cgrid[c],fc_ij)*t;
+ }
+ }
+ }
+}
+
+size_t DM_EAD::gen_atomic_orbitals(int Lmax) {
+ orbitals.resize(Lmax+1);
+ fac.resize(Lmax+1);
+
+ size_t count = 0;
+ for (int L = 0; L <= Lmax; ++L) {
+ for (int lx = 0; lx <= L; ++lx) {
+ for (int ly = 0; ly <= L; ++ly) {
+ for (int lz = 0; lz <= L; ++lz) {
+ if (lx+ly+lz == L) {
+ std::vector<int> o = {lx, ly, lz};
+ orbitals[L].push_back(o);
+ const double f = fact(L)/(fact(lx)*fact(ly)*fact(lz))/my_pow(4*rc,L);
+ fac[L].push_back(f);
+ count++;
+ }
+ }
+ }
+ }
+ }
+ return count;
+}
+double DM_EAD::fact(int n)
+{
+ double f=1.0;
+ while (n) {
+ f *= n;
+ --n;
+ }
+ return f;
+}
+double DM_EAD::my_pow(double x, int n){
+ double r = 1.0;
+
+ while(n > 0){
+ r *= x;
+ --n;
+ }
+ return r;
+}
+size_t DM_EAD::rhoi_size() {
+ return rhoisize;
+}
+size_t DM_EAD::rhoip_size() {
+ return rhoisize;
+}
diff --git a/descriptors/dm/dm_ead.h b/descriptors/dm/dm_ead.h
new file mode 100644
index 0000000000000000000000000000000000000000..3c85158f6168401f243715f3f6c0ac7f22855aca
--- /dev/null
+++ b/descriptors/dm/dm_ead.h
@@ -0,0 +1,89 @@
+#ifndef DM_EAD_H
+#define DM_EAD_H
+#include "dm_base.h"
+
+/** Embedded Atom Descriptor
+ *
+ * \f[
+ * V_i^{L,\eta,r_s} = \sum_{l_x,l_y,l_z}^{l_x+l_y+l_z=L} \frac{L!}{l_x!l_y!l_z!}
+ * \Big( \rho_i^{\eta,r_s,l_x,l_y,l_z} \Big)^2
+ * \f]
+ *
+ * where density \f$ \rho \f$ is calculated using Gaussian Type Orbitals:
+ * \f[
+ * \rho_i^{\eta,r_s,l_x,l_y,l_z} = \sum_{j \neq i} x_{ij}^{l_x}y_{ij}^{l_y}z_{ij}^{l_z}
+ * \exp{\Big(-\eta(r_{ij}-r_s)^2\Big)}f_c(r_{ij})
+ * \f]
+ *
+ * \ref CGRIDMB parameters control position \f$ r_s \f$ of the gaussian basis function.
+ *
+ * \ref SGRIDMB parameters control width \f$ \eta \f$ of the gaussian basis function.
+ *
+ * \ref AGRIDMB parameter specify maximum value for the angular momentum \f$L_{max}\f$.
+ *
+ * e.g. \f$L_{max}=2\f$ will calculate descriptors with \f$ L=0,1,2 \f$ (s,p,d orbitals).
+ *
+ * More information about this descriptor:
+ *
+ * <div class="csl-entry">Zhang, Y., Hu, C.,Jiang, B.
+ * (2019). Embedded atom neural network potentials: efficient and accurate
+ * machine learning with a physically inspired representation.
+ * <i>Journal of Physical Chemistry Letters</i>, <i>10</i>(17),
+ * 4962–4967. https://doi.org/10.1021/acs.jpclett.9b02037</div>
+ *
+ * Required Config keys:
+ * \ref INITMB \ref CGRIDMB \ref SGRIDMB \ref AGRIDMB
+ */
+class DM_EAD: public DM_Base {
+ private:
+ int Lmax=0;
+ size_t rhoisize=0;
+ size_t s=0;
+ std::string lab="DM_EAD";
+ v_type sgrid;
+ v_type cgrid;
+ std::vector<std::vector<std::vector<int>>> orbitals;
+ std::vector<std::vector<double>> fac;
+ size_t gen_atomic_orbitals(int L);
+ double fact(int n);
+ double my_pow(double,int);
+ int verbose;
+ double rc;
+
+ public:
+ DM_EAD(Config&);
+ void calc_aed(
+ rho_rtype rho,
+ aed_rtype aed);
+ int calc_dXijdri_dXjidri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj);
+ int calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+ void init_rhoi(rho_rtype rhoi);
+ void calc_rho(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rho);
+ size_t rhoi_size();
+ size_t rhoip_size();
+};
+#endif
diff --git a/descriptors/dm/dm_eam.cpp b/descriptors/dm/dm_eam.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..56c2b8533b059d72230fc9d32dbfb49f3c97a4ce
--- /dev/null
+++ b/descriptors/dm/dm_eam.cpp
@@ -0,0 +1,222 @@
+#include "dm_eam.h"
+
+Registry<DM_Base,Config&>::Register<DM_EAM> DM_EAM_1( "DM_EAM" );
+
+DM_EAM::DM_EAM(Config &c):
+ verbose(c.get<int>("VERBOSE"))
+{
+ if (!c.get<bool>("INITMB")) {
+ s=0;
+ return;
+ }
+ ef.file_path = c("SETFL")[0];
+ read_setfl();
+
+ if (abs(ef.rcut - c.get<double>("RCUTMB")) > std::numeric_limits<double>::min()) {
+ if (verbose) {
+ std::cout << "Config file cutoff and setfl file cutoff differ: "
+ << c.get<double>("RCUTMB") << " " << ef.rcut << std::endl;
+ std::cout << "Enforcing SETFL file cutoff: " << ef.rcut << std::endl;
+ }
+ c.remove("RCUTMB");
+ c.add("RCUTMB", ef.rcut);
+ c.postprocess();
+ }
+
+
+ frho_spline.resize(ef.nrho+1, std::vector<double>(7));
+ rhor_spline.resize(ef.nr+1, std::vector<double>(7));
+ z2r_spline.resize(ef.nr+1, std::vector<double>(7));
+
+
+ gen_splines(ef.nrho, ef.drho, ef.frho, frho_spline);
+ gen_splines(ef.nr, ef.dr, ef.rhor, rhor_spline);
+ gen_splines(ef.nr, ef.dr, ef.z2r, z2r_spline);
+
+}
+
+void DM_EAM::calc_aed(
+ rho_rtype rho,
+ aed_rtype aed)
+{
+ //std::cout << "rho[i]: " << rho(0) << std::endl;
+ double p = rho(0)*ef.rdrho + 1.0;
+ int m = static_cast<int> (p);
+ m = std::max(1,std::min(m,ef.nrho-1));
+ p -= m;
+ p = std::min(p,1.0);
+ //std::vector<double> coeff = frho_spline[m];
+ double phi = ((frho_spline[m][3]*p + frho_spline[m][4])*p + frho_spline[m][5])*p + frho_spline[m][6];
+
+ rho(1) = (frho_spline[m][0]*p + frho_spline[m][1])*p + frho_spline[m][2]; // lammps fp[i]
+ //std::cout << "fp[i]: " << rho(1) << std::endl;
+
+ if (rho(0) > ef.rhomax) phi += rho(1) * (rho(0)-ef.rhomax);
+ //phi *= scale[type[i]][type[i]];
+
+ aed(fidx) += phi;
+
+
+}
+int DM_EAM::calc_dXijdri_dXjidri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &,
+ const double,
+ const double,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double,
+ const double)
+{
+ double r = rij;
+ //double recip = 1.0/r;
+ double p = r*ef.rdr + 1.0;
+ int m = static_cast<int> (p);
+ m = std::min(m,ef.nr-1);
+ p -= m;
+ p = std::min(p,1.0);
+ //std::vector<double> coeff = rhor_spline[m];
+ double rhoip = (rhor_spline[m][0]*p + rhor_spline[m][1])*p + rhor_spline[m][2];
+ double rhojp = (rhor_spline[m][0]*p + rhor_spline[m][1])*p + rhor_spline[m][2];
+
+ // here rho_i[1] is a derivative of rho_i at local atom, similarly for rho_j
+ //double psip = rho_i[1]*rhojp + rho_j[1]*rhoip;
+ double psip = rhoi(1)*rhojp + rhoj(1)*rhoip;
+
+ fd_ij(fidx,0) = psip; // requires *delij*recip in the main loop
+ return 0;
+
+}
+int DM_EAM::calc_dXijdri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &,
+ const double,
+ const double,
+ rho_rtype rhoi,
+ fd_type &fd_ij)
+{
+ double r = rij;
+ //double recip = 1.0/r;
+ double p = r*ef.rdr + 1.0;
+ int m = static_cast<int> (p);
+ m = std::min(m,ef.nr-1);
+ p -= m;
+ p = std::min(p,1.0);
+ double rhojp = (rhor_spline[m][0]*p + rhor_spline[m][1])*p + rhor_spline[m][2];
+
+ // here rho_i[1] is a derivative of rho_i at local atom
+ double psip = rhoi(1)*rhojp;
+
+ fd_ij(fidx,0) = psip; // requires *delij*recip in the main loop
+ return 0;
+
+}
+size_t DM_EAM::size() {
+ return s;
+}
+std::string DM_EAM::label() {
+ return lab;
+}
+
+void DM_EAM::read_setfl()
+{
+ std::string line;
+ std::ifstream in_file(ef.file_path);
+ if (!in_file.good())
+ throw std::runtime_error("SETFL file does not exists.\n");
+
+ if (in_file.is_open()) {
+ if (verbose)
+ std::cout << "<DM_EAM> Reading setfl: " << ef.file_path << std::endl;
+ // skip ficgridt 3 comment lines
+ getline (in_file,line);
+ getline (in_file,line);
+ getline (in_file,line);
+ // skip number of types and types
+ getline (in_file,line);
+ // read 5th line
+ in_file >> ef.nrho >> ef.drho >> ef.nr >> ef.dr >> ef.rcut;
+ in_file >> ef.atomic_number >> ef.atomic_mass >> ef.lattice_param >> ef.lattice;
+ ef.rdr = 1.0/ef.dr;
+ ef.rdrho = 1.0/ef.drho;
+ // prepare arrays
+ ef.frho.resize(ef.nrho);
+ ef.rhor.resize(ef.nr);
+ ef.z2r.resize(ef.nr);
+ // read all data
+ for (int i=0; i<ef.nrho; ++i) in_file >> ef.frho[i];
+ for (int i=0; i<ef.nr; ++i) in_file >> ef.rhor[i];
+ for (int i=0; i<ef.nr; ++i) in_file >> ef.z2r[i];
+ in_file.close();
+ }
+ else {
+ std::cout << "<DM_EAM> Unable to open file: " << ef.file_path << std::endl;;
+ }
+ // get max value of rho
+ //ef.rhomax = *std::max_element(ef.rhor.begin(), ef.rhor.end());
+ ef.rhomax = (ef.nrho-1) * ef.drho;
+
+}
+void DM_EAM::gen_splines(int &n, double &delta, std::vector<double> &f,
+ std::vector<std::vector<double>> &spline)
+{
+ // in lammps f is n+1, here is size n
+ for (int m=1; m<=n; m++) spline[m][6] = f[m-1];
+
+ spline[1][5] = spline[2][6] - spline[1][6];
+ spline[2][5] = 0.5 * (spline[3][6]-spline[1][6]);
+ spline[n-1][5] = 0.5 * (spline[n][6]-spline[n-2][6]);
+ spline[n][5] = spline[n][6] - spline[n-1][6];
+
+ for (int m = 3; m <= n-2; m++)
+ spline[m][5] = ((spline[m-2][6]-spline[m+2][6]) +
+ 8.0*(spline[m+1][6]-spline[m-1][6])) / 12.0;
+
+ for (int m = 1; m <= n-1; m++) {
+ spline[m][4] = 3.0*(spline[m+1][6]-spline[m][6]) - 2.0*spline[m][5] - spline[m+1][5];
+ spline[m][3] = spline[m][5] + spline[m+1][5] - 2.0*(spline[m+1][6]-spline[m][6]);
+ }
+
+ spline[n][4] = 0.0;
+ spline[n][3] = 0.0;
+
+ for (int m = 1; m <= n; m++) {
+ spline[m][2] = spline[m][5]/delta;
+ spline[m][1] = 2.0*spline[m][4]/delta;
+ spline[m][0] = 3.0*spline[m][3]/delta;
+ }
+}
+void DM_EAM::init_rhoi(rho_rtype rhoi)
+{
+ rhoi.resize(2);
+ rhoi.setZero();
+ //rho.resize(1);
+ //rho[0] = Eigen::Matrix<double,1,2>(1,2); // [0] - rho, [1] - rho prime aka fp[i]
+}
+void DM_EAM::calc_rho(
+ const double rij,
+ const double,
+ const double,
+ const Eigen::Vector3d &,
+ rho_rtype rho)
+{
+ double r = rij;
+ double p = r*ef.rdr + 1.0;
+ int m = static_cast<int> (p);
+ m = std::min(m,ef.nr-1);
+ p -= m;
+ p = std::min(p,1.0);
+ std::vector<double> coeff = rhor_spline[m];
+ // rho[i]
+ rho(0) += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+
+}
+size_t DM_EAM::rhoi_size() {
+ return 1;
+}
+size_t DM_EAM::rhoip_size() {
+ return 1;
+}
diff --git a/descriptors/dm/dm_eam.h b/descriptors/dm/dm_eam.h
new file mode 100644
index 0000000000000000000000000000000000000000..e2808544f128baa6d311f84041bddace1438bb9e
--- /dev/null
+++ b/descriptors/dm/dm_eam.h
@@ -0,0 +1,93 @@
+#ifndef DM_EAM_H
+#define DM_EAM_H
+#include "dm_base.h"
+
+/** \brief many-body part for the Embedded Atom Method descriptor.
+ *
+ * \f[
+ * V_i = F\Bigg(\sum_{j \neq i} \rho(r_{ij}) \Bigg)
+ * \f]
+ *
+ * This descriptor will load tabulated values for the density \f$ \rho \f$
+ * and embedded energy \f$ F \f$ from the provided \ref SETFL file.
+ *
+ * This descriptor is usually used together with the two-body descriptor \ref D2_EAM
+ * although this is not required and user can mix it with any other descriptors
+ * or use it on its own.
+ *
+ * This descriptor will enforce cutoff distance as specified in a \ref SETFL file.
+ * Set \ref RCUTMB to the same value to suppress the warning message.
+ *
+ * Required Config keys:
+ * \ref INITMB \ref SETFL
+ */
+class DM_EAM: public DM_Base {
+ private:
+ int verbose;
+ struct eam_file {
+ std::string file_path;
+ std::vector<double> frho;
+ std::vector<double> rhor;
+ std::vector<double> z2r;
+ int nrho=0;
+ double drho=0;
+ int nr;
+ double dr;
+ double rdr;
+ double rdrho;
+ double rcut;
+ int atomic_number;
+ double atomic_mass;
+ double lattice_param;
+ std::string lattice;
+ double rhomax;
+
+ };
+ eam_file ef;
+
+ std::vector<std::vector<double>> frho_spline;
+ std::vector<std::vector<double>> rhor_spline;
+ std::vector<std::vector<double>> z2r_spline;
+
+ size_t s=1;
+ std::string lab="DM_EAM";
+ void read_setfl();
+ void gen_splines(int &n, double &delta, std::vector<double> &f, std::vector<std::vector<double>> &spline);
+
+ public:
+ DM_EAM(Config&);
+ void calc_aed(
+ rho_rtype rho,
+ aed_rtype aed);
+ int calc_dXijdri_dXjidri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj);
+ int calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+ void init_rhoi(rho_rtype rhoi);
+ void calc_rho(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rho);
+ size_t rhoi_size();
+ size_t rhoip_size();
+};
+#endif
diff --git a/descriptors/dm/dm_mEAD.cpp b/descriptors/dm/dm_mEAD.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..dd2a503ba93d6e63fa774f45c8a1baee4144c5b9
--- /dev/null
+++ b/descriptors/dm/dm_mEAD.cpp
@@ -0,0 +1,318 @@
+#include "dm_mEAD.h"
+#include "../d_basis_functions.h"
+#include <vector>
+#include <iomanip>
+
+
+Registry<DM_Base,Config&>::Register<DM_mEAD<F_RLR>> DM_mEAD_1( "DM_mEAD1" );
+
+/* Index ii iterates rhoi array
+ * rhoi array is split 50/50 into density and derivative
+ * of the embedding function
+ */
+template <typename F>
+DM_mEAD<F>::DM_mEAD(Config &config):
+ verbose(config.get<int>("VERBOSE")),
+ rc(config.get<double>("RCUTMB")),
+ emb(config)
+{
+
+ if (!config.get<bool>("INITMB")) return;
+
+ get_grid(config,"CGRIDMB",cgrid);
+ get_grid(config,"SGRIDMB",sgrid);
+
+ if (cgrid.size()!=sgrid.size()) {
+ throw std::runtime_error("SGRID2B and CGRID2B arrays differ in size.\n");
+ }
+
+ Lmax = config.get<int>("AGRIDMB");
+ s=sgrid.size()*(Lmax+1);
+ rhoisize = gen_atomic_orbitals(Lmax);
+ rhoisize *= sgrid.size();
+
+ if (verbose) {
+ std::cout << std::endl;
+ std::cout << "SGRID (SGRIDMB): " << sgrid.size() << std::endl;
+ for (auto e:sgrid) std::cout << e << " ";
+ std::cout << std::endl;
+
+ std::cout << "CGRID (CGRIDMB): " << cgrid.size() << std::endl;
+ for (auto L:cgrid) std::cout << L << " ";
+ std::cout << std::endl;
+ std::cout << "rhoisize: " << rhoisize << std::endl;
+ }
+}
+
+template <typename F>
+void DM_mEAD<F>::calc_aed(
+ rho_rtype rhoi,
+ aed_rtype aed)
+{
+ size_t Lshift = 0;
+ size_t ii=0;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const double f = fac[L][o];
+ for (size_t c=0; c<cgrid.size(); c++) {
+ aed(c+Lshift+fidx) += f*emb.calc_F(rhoi(ii),c);
+ rhoi(ii+rhoisize) = f*emb.calc_dF(rhoi(ii),c);
+ ii++;
+ }
+ }
+ Lshift += sgrid.size();
+ }
+}
+template <typename F>
+int DM_mEAD<F>::calc_dXijdri_dXjidri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj)
+{
+ const double x = vec_ij[0];
+ const double y = vec_ij[1];
+ const double z = vec_ij[2];
+
+ size_t Lshift = 0;
+
+ size_t ii=rhoisize;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const int lx = orbitals[L][o][0];
+ const int ly = orbitals[L][o][1];
+ const int lz = orbitals[L][o][2];
+
+ const double powxlxij = my_pow(x,lx);
+ const double powylyij = my_pow(y,ly);
+ const double powzlzij = my_pow(z,lz);
+
+ const double powxlxji = my_pow(-x,lx);
+ const double powylyji = my_pow(-y,ly);
+ const double powzlzji = my_pow(-z,lz);
+
+ double txij=0.0,tyij=0.0,tzij=0.0;
+ double txji=0.0,tyji=0.0,tzji=0.0;
+ if (lx!=0) {
+ txij = lx*my_pow(x,lx-1)*powylyij*powzlzij;
+ txji = lx*my_pow(-x,lx-1)*powylyji*powzlzji;
+ }
+ if (ly!=0) {
+ tyij = ly*my_pow(y,ly-1)*powxlxij*powzlzij;
+ tyji = ly*my_pow(-y,ly-1)*powxlxji*powzlzji;
+ }
+ if (lz!=0) {
+ tzij = lz*my_pow(z,lz-1)*powylyij*powxlxij;
+ tzji = lz*my_pow(-z,lz-1)*powylyji*powxlxji;
+ }
+ const double pqrij = powxlxij*powylyij*powzlzij;
+ const double pqrji = powxlxji*powylyji*powzlzji;
+
+ for (size_t c=0; c<cgrid.size(); c++) {
+ double bf= G(rij,sgrid[c],cgrid[c],fc_ij);
+ double dbf= dG(rij,sgrid[c],cgrid[c],fc_ij,fcp_ij);
+
+ const double term2ijx=pqrij*dbf*x/rij;
+ const double term2ijy=pqrij*dbf*y/rij;
+ const double term2ijz=pqrij*dbf*z/rij;
+
+ const double term2jix=-pqrji*dbf*x/rij;
+ const double term2jiy=-pqrji*dbf*y/rij;
+ const double term2jiz=-pqrji*dbf*z/rij;
+
+ const double term1ijx = bf*txij;
+ const double term1ijy = bf*tyij;
+ const double term1ijz = bf*tzij;
+
+ const double term1jix = bf*txji;
+ const double term1jiy = bf*tyji;
+ const double term1jiz = bf*tzji;
+
+ fd_ij(c+Lshift+fidx,0) +=
+ rhoi(ii)*(term1ijx + term2ijx)*wj
+ -rhoj(ii)*(term1jix + term2jix)*wi
+ ;
+
+ fd_ij(c+Lshift+fidx,1) +=
+ rhoi(ii)*(term1ijy + term2ijy)*wj
+ -rhoj(ii)*(term1jiy + term2jiy)*wi
+ ;
+
+ fd_ij(c+Lshift+fidx,2) +=
+ rhoi(ii)*(term1ijz + term2ijz)*wj
+ -rhoj(ii)*(term1jiz + term2jiz)*wi
+ ;
+ ii++;
+ }
+ }
+ Lshift+=sgrid.size();
+ }
+
+ return 1;
+
+}
+template <typename F>
+int DM_mEAD<F>::calc_dXijdri(
+ const double rij,
+ const double,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij)
+{
+ const double x = vec_ij[0];
+ const double y = vec_ij[1];
+ const double z = vec_ij[2];
+
+ size_t Lshift = 0;
+
+ size_t ii=rhoisize;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const int lx = orbitals[L][o][0];
+ const int ly = orbitals[L][o][1];
+ const int lz = orbitals[L][o][2];
+
+ const double powxlxij = my_pow(x,lx);
+ const double powylyij = my_pow(y,ly);
+ const double powzlzij = my_pow(z,lz);
+
+ double txij=0.0,tyij=0.0,tzij=0.0;
+ if (lx!=0) {
+ txij = lx*my_pow(x,lx-1)*powylyij*powzlzij;
+ }
+ if (ly!=0) {
+ tyij = ly*my_pow(y,ly-1)*powxlxij*powzlzij;
+ }
+ if (lz!=0) {
+ tzij = lz*my_pow(z,lz-1)*powylyij*powxlxij;
+ }
+ const double pqrij = powxlxij*powylyij*powzlzij;
+
+ for (size_t c=0; c<cgrid.size(); c++) {
+ double bf= G(rij,sgrid[c],cgrid[c],fc_ij);
+ double dbf= dG(rij,sgrid[c],cgrid[c],fc_ij,fcp_ij);
+
+ const double term2ijx=pqrij*dbf*x/rij;
+ const double term2ijy=pqrij*dbf*y/rij;
+ const double term2ijz=pqrij*dbf*z/rij;
+
+ const double term1ijx = bf*txij;
+ const double term1ijy = bf*tyij;
+ const double term1ijz = bf*tzij;
+
+ fd_ij(c+Lshift+fidx,0) +=
+ rhoi(ii)*(term1ijx + term2ijx)
+ ;
+
+ fd_ij(c+Lshift+fidx,1) +=
+ rhoi(ii)*(term1ijy + term2ijy)
+ ;
+
+ fd_ij(c+Lshift+fidx,2) +=
+ rhoi(ii)*(term1ijz + term2ijz)
+ ;
+ ii++;
+ }
+ }
+ Lshift+=sgrid.size();
+ }
+
+ return 1;
+
+}
+template <typename F>
+size_t DM_mEAD<F>::size() {
+ return s;
+}
+template <typename F>
+std::string DM_mEAD<F>::label() {
+ return lab;
+}
+
+template <typename F>
+void DM_mEAD<F>::init_rhoi(rho_rtype rhoi)
+{
+ rhoi.resize(2*rhoisize);
+ rhoi.setZero();
+}
+template <typename F>
+void DM_mEAD<F>::calc_rho(
+ const double rij,
+ const double,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rhoi)
+{
+ size_t ii=0;
+ for (int L=0; L<=Lmax; ++L) {
+ for (size_t o=0; o<orbitals[L].size(); ++o) {
+ const size_t lx = orbitals[L][o][0];
+ const size_t ly = orbitals[L][o][1];
+ const size_t lz = orbitals[L][o][2];
+ double t = (my_pow(vec_ij[0],lx)*my_pow(vec_ij[1],ly)*my_pow(vec_ij[2],lz));
+ for (size_t c=0; c<cgrid.size(); c++) {
+ rhoi(ii++) += G(rij,sgrid[c],cgrid[c],fc_ij)*t;
+ }
+ }
+ }
+}
+
+template <typename F>
+size_t DM_mEAD<F>::gen_atomic_orbitals(int Lmax) {
+ orbitals.resize(Lmax+1);
+ fac.resize(Lmax+1);
+
+ size_t count = 0;
+ for (int L = 0; L <= Lmax; ++L) {
+ for (int lx = 0; lx <= L; ++lx) {
+ for (int ly = 0; ly <= L; ++ly) {
+ for (int lz = 0; lz <= L; ++lz) {
+ if (lx+ly+lz == L) {
+ std::vector<int> o = {lx, ly, lz};
+ orbitals[L].push_back(o);
+ const double f = fact(L)/(fact(lx)*fact(ly)*fact(lz))/my_pow(4*rc,L);
+ fac[L].push_back(f);
+ count++;
+ }
+ }
+ }
+ }
+ }
+ return count;
+}
+template <typename F>
+double DM_mEAD<F>::fact(int n)
+{
+ double f=1.0;
+ while (n) {
+ f *= n;
+ --n;
+ }
+ return f;
+}
+template <typename F>
+double DM_mEAD<F>::my_pow(double x, int n){
+ double r = 1.0;
+
+ while(n > 0){
+ r *= x;
+ --n;
+ }
+ return r;
+}
+template <typename F>
+size_t DM_mEAD<F>::rhoi_size() {
+ return rhoisize;
+}
+template <typename F>
+size_t DM_mEAD<F>::rhoip_size() {
+ return rhoisize;
+}
diff --git a/descriptors/dm/dm_mEAD.h b/descriptors/dm/dm_mEAD.h
new file mode 100644
index 0000000000000000000000000000000000000000..27f859c600c9fa4c16bc8cf725cd9426a1beb4fc
--- /dev/null
+++ b/descriptors/dm/dm_mEAD.h
@@ -0,0 +1,146 @@
+#ifndef DM_mEAD_H
+#define DM_mEAD_H
+#include "dm_base.h"
+class F_Base {
+ public:
+ virtual double calc_F(double rho, size_t c)=0;
+ virtual double calc_dF(double rho, size_t c)=0;
+};
+
+/** Embedding function of the form: sx*log(cx)
+ *
+ * REQUIRED KEYS:
+ *
+ * SEMBFUNC, CEMBFUNC
+ * control depth and x-intercept of an embedding function respectively.
+ * The number of keys must match those in the mEAD descriptor.
+ *
+ */
+class F_RLR: public F_Base {
+ private:
+ Config *config=nullptr;
+ v_type sgrid; // controls depth
+ v_type cgrid; // control x-intercep at 1/c
+ public:
+ F_RLR(Config &conf):
+ config(&conf)
+ {
+ if (!config->get<bool>("INITMB")) return;
+ D_Base::get_grid(*config,"SEMBFUNC",sgrid);
+ D_Base::get_grid(*config,"CEMBFUNC",cgrid);
+
+ }
+ F_RLR():
+ config(nullptr)
+ {
+ }
+ double calc_F(double rho,size_t c) {
+ if (rho>0)
+ return sgrid[c]*rho*log(cgrid[c]*rho);
+ return 0;
+ };
+ double calc_dF(double rho,size_t c) {
+ if (rho>0)
+ return sgrid[c]*(log(cgrid[c]*rho)+1);
+ return sgrid[c];
+ };
+ //~F_RLR() {
+ // Do not delete config!
+ //}
+};
+
+/** Modified Embedded Atom Descriptor
+ *
+ * REQUIRED KEYS: SCGRIDMB, CGRIDMB
+ * + REQUIRED KEYS BY THE EMBEDDING FUNCTION
+ *
+ *
+ * TODO Below description is for EAD not mEAD
+ * This descriptor has mathematical form very similar to EAD.
+ * It allow: for different embedding functions to be used.
+ *
+ * \f[
+ * V_i^{L,\eta,r_s} = \sum_{l_x,l_y,l_z}^{l_x+l_y+l_z=L} \frac{L!}{l_x!l_y!l_z!}
+ * \Big( \rho_i^{\eta,r_s,l_x,l_y,l_z} \Big)^2
+ * \f]
+ *
+ * where density \f$ \rho \f$ is calculated using Gaussian Type Orbitals:
+ * \f[
+ * \rho_i^{\eta,r_s,l_x,l_y,l_z} = \sum_{j \neq i} x_{ij}^{l_x}y_{ij}^{l_y}z_{ij}^{l_z}
+ * \exp{\Big(-\eta(r_{ij}-r_s)^2\Big)}f_c(r_{ij})
+ * \f]
+ *
+ * \ref CGRIDMB parameters control position \f$ r_s \f$ of the gaussian basis function.
+ *
+ * \ref SGRIDMB parameters control width \f$ \eta \f$ of the gaussian basis function.
+ *
+ * \ref AGRIDMB parameter specify maximum value for the angular momentum \f$L_{max}\f$.
+ *
+ * e.g. \f$L_{max}=2\f$ will calculate descriptors with \f$ L=0,1,2 \f$ (s,p,d orbitals).
+ *
+ * More information about this descriptor:
+ *
+ * <div class="csl-entry">Zhang, Y., Hu, C.,Jiang, B.
+ * (2019). Embedded atom neural network potentials: efficient and accurate
+ * machine learning with a physically inspired representation.
+ * <i>Journal of Physical Chemistry Letters</i>, <i>10</i>(17),
+ * 4962–4967. https://doi.org/10.1021/acs.jpclett.9b02037</div>
+ *
+ * Required Config keys:
+ * \ref INITMB \ref CGRIDMB \ref SGRIDMB \ref AGRIDMB
+ */
+template <typename F>
+class DM_mEAD: public DM_Base, public F {
+ private:
+ int Lmax=0;
+ size_t rhoisize=0;
+ size_t s=0;
+ std::string lab="DM_mEAD";
+ v_type sgrid;
+ v_type cgrid;
+ std::vector<std::vector<std::vector<int>>> orbitals;
+ std::vector<std::vector<double>> fac;
+ size_t gen_atomic_orbitals(int L);
+ double fact(int n);
+ double my_pow(double,int);
+ int verbose;
+ double rc;
+ F emb;
+
+ public:
+ DM_mEAD(Config&);
+ void calc_aed(
+ rho_rtype rho,
+ aed_rtype aed);
+ int calc_dXijdri_dXjidri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ rho_rtype rhoj,
+ fd_type &fd_ij,
+ const double wi,
+ const double wj);
+ int calc_dXijdri(
+ const double rij,
+ const double rij_sq,
+ const Eigen::Vector3d &vec_ij,
+ const double fc_ij,
+ const double fcp_ij,
+ rho_rtype rhoi,
+ fd_type &fd_ij);
+ size_t size();
+ std::string label();
+ void init_rhoi(rho_rtype rhoi);
+ void calc_rho(
+ const double rij,
+ const double rij_sq,
+ const double fc_ij,
+ const Eigen::Vector3d &vec_ij,
+ rho_rtype rho);
+ size_t rhoi_size();
+ size_t rhoip_size();
+};
+#endif
diff --git a/ea.h b/ea.h
new file mode 100644
index 0000000000000000000000000000000000000000..528fe973a7a9a862469ed32930cdd26fe2c8f3c6
--- /dev/null
+++ b/ea.h
@@ -0,0 +1,90 @@
+#ifndef M_EA_H
+#define M_EA_H
+
+#include "../LIBS/Eigen/Dense"
+#include "../CORE/config/config.h"
+#include "../CORE/typedefs.h"
+
+class EA {
+ private:
+ Config &config;
+ int verbose;
+ public:
+ EA(Config &c):
+ config(c),
+ verbose(c.get<int>("VERBOSE"))
+ {}
+
+ int run (const mat &Phi,const vec &T,
+ const mat &S, const mat &U,const mat &V,
+ double &alpha, double &beta) {
+
+ // Phi = USV^T
+ size_t N=Phi.rows();
+ size_t maxsteps=40;
+ const double EPS2 = 5e-12;
+ int convergence_status = 1;
+ double test1 = 2.0*EPS2;
+ double test2 = 2.0*EPS2;
+ double alpha_old = alpha;
+ double beta_old = beta;
+ size_t counter=0;
+ double lambda = alpha/beta;
+
+ // square of the S matrix contains all eigenvalues
+ // of Phi^T*Phi = (VSU^T)(USV^T) = VS^2V^T = S^2
+ const vec eval = S.array().pow(2);
+
+
+ if (verbose) printf("%9s %9s %9s %9s %9s %9s\n","del alpha", "del beta", "alpha", "beta", "lambda", "gamma");
+ // D is a filter factors matrix
+ mat D(S.rows(),S.rows());
+ D.setZero();
+ int outprec=config.get<int>("OUTPREC");
+
+ while (test1>EPS2 || test2>EPS2) {
+
+ for (long int i=0; i<S.rows(); ++i) {
+ double s = S(i,0);
+ if (s>std::numeric_limits<double>::min())
+ D(i,i) = (s)/(s*s+lambda);
+ else
+ D(i,i) = 0.0;
+ }
+ // regularised least square problem (Tikhonov Regularization):
+ vec m_n = V*D*U.transpose()*T;
+
+ double gamma = 0.0;
+ for (long int j=0; j<S.rows(); j++) {
+ // exlude eigenvalues which are smaller than ~1e-10
+ // for numerical stability b/c matrix is not regularized
+ if (beta*eval(j) > std::numeric_limits<double>::min()) {
+ gamma += (beta*eval(j))/(beta*eval(j)+alpha);
+ }
+ }
+ alpha = gamma / m_n.dot(m_n);
+
+ double temp = (T-Phi*m_n).dot(T-Phi*m_n);
+ beta = (static_cast<double>(N)-gamma)/temp;
+
+ test1 = fabs(alpha_old-alpha)/alpha;
+ test2 = fabs(beta_old-beta)/beta;
+ if (verbose) printf("%.3e %.3e %.3e %.3e %.3e %.3e\n",test1, test2,alpha, beta, alpha/beta, gamma);
+
+ lambda = alpha/beta;
+ alpha_old = alpha;
+ beta_old = beta;
+
+ if (++counter>maxsteps) {
+ convergence_status = 0;
+ break;
+ }
+ }
+ if (convergence_status && verbose) std::cout << "Number of steps for convergence: " + to_string(counter,outprec) << std::endl;
+ if (!convergence_status && verbose) std::cout << "Max number of steps reached: " + to_string(counter,outprec) << std::endl;
+ return convergence_status;
+
+ }
+};
+#endif
+
diff --git a/ekm.h b/ekm.h
new file mode 100644
index 0000000000000000000000000000000000000000..5f60bc487b52622259d5fc2a77e7a592ccdb4c35
--- /dev/null
+++ b/ekm.h
@@ -0,0 +1,50 @@
+#ifndef EKM_H
+#define EKM_H
+
+#include "../LIBS/Eigen/Dense"
+#include "../CORE/typedefs.h"
+
+template <typename K>
+class EKM {
+ private:
+ K kernel;
+ public:
+ Eigen::MatrixXd KK;
+ template <typename T>
+ EKM(T &t):
+ kernel(t)
+ {}
+ void configure(aeds_type &basis) {
+ // KK is temp matrix
+ KK = Eigen::MatrixXd(basis.size(), basis.size());
+
+ for (long int i=0; i<KK.rows(); i++) {
+ for (long int j=0; j<KK.cols(); j++) {
+ double val = kernel(basis.col(i),basis.col(j));
+ KK(i,j) = val;
+ }
+ }
+
+ // KK is a kernel matrix
+ // we compute cholesky for KK^-1
+ // then KK = L^-1
+ KK = KK.inverse();
+ Eigen::LLT<Eigen::MatrixXd> llt(KK);
+ KK = llt.matrixL();
+ KK.transposeInPlace();
+ }
+
+ void project(phi_type &Phi) {
+ for (long int i=0; i<Phi.rows(); ++i) {
+ Phi.row(i)= KK*Phi.row(i).transpose();
+ }
+ }
+ aed_type project(aed_rtype aed, aeds_type &basis) {
+ aed_type temp(basis.size());
+ for (long int b=0; b<basis.cols(); b++) {
+ temp(b) = kernel(aed,basis[b]);
+ }
+ return KK*temp;
+ }
+};
+#endif
diff --git a/functions/basis_functions/bf_all.h b/functions/basis_functions/bf_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..6cef8e3498dc1527d5091e2cc9f4ab3eddcceb45
--- /dev/null
+++ b/functions/basis_functions/bf_all.h
@@ -0,0 +1,2 @@
+#include "bf_linear.h"
+#include "bf_polynomial2.h"
diff --git a/functions/basis_functions/bf_base.cpp b/functions/basis_functions/bf_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..82eb1340d45bfd644bfdc9ce1291dffd12c74a1f
--- /dev/null
+++ b/functions/basis_functions/bf_base.cpp
@@ -0,0 +1,2 @@
+#include "bf_base.h"
+BF_Base::~BF_Base() {}
diff --git a/functions/basis_functions/bf_base.h b/functions/basis_functions/bf_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..d76cd7194c918aaccedf85cf9f6ef1e8c3ae9384
--- /dev/null
+++ b/functions/basis_functions/bf_base.h
@@ -0,0 +1,18 @@
+#ifndef BASIS_FUNCTIONS_H
+#define BASIS_FUNCTIONS_H
+
+#include "../../../CORE/typedefs.h"
+#include "../function_base.h"
+
+#include <iostream>
+
+/** \brief Base class to be used by all basis functions */
+struct BF_Base: public virtual Function_Base {
+ virtual ~BF_Base();
+ virtual double epredict(const t_type &, const aed_rctype& )=0;
+ virtual double fpredict(const t_type &, const fd_type &,
+ const aed_rctype& , const size_t )=0;
+ virtual force_type fpredict(const t_type &, const fd_type &,
+ const aed_rctype& )=0;
+};
+#endif
diff --git a/functions/basis_functions/bf_linear.cpp b/functions/basis_functions/bf_linear.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..714bab185c773131fdc41233e8c7c64fd5430fd3
--- /dev/null
+++ b/functions/basis_functions/bf_linear.cpp
@@ -0,0 +1,25 @@
+#include "bf_linear.h"
+
+Registry<Function_Base>::Register<BF_Linear> BF_Linear_1( "BF_Linear" );
+Registry<Function_Base,Config&>::Register<BF_Linear> BF_Linear_2( "BF_Linear" );
+
+
+BF_Linear::BF_Linear() {}
+BF_Linear::BF_Linear(const Config &c): Function_Base(c) {}
+std::string BF_Linear::get_label() {
+ return label;
+}
+double BF_Linear::epredict(const t_type &weights, const aed_rctype& aed)
+{
+ return weights.transpose()*aed;
+}
+double BF_Linear::fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& , const size_t k)
+{
+ return -fdij.col(k).transpose() * weights;
+}
+force_type BF_Linear::fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& )
+{
+ return -fdij.transpose() * weights;
+}
diff --git a/functions/basis_functions/bf_linear.h b/functions/basis_functions/bf_linear.h
new file mode 100644
index 0000000000000000000000000000000000000000..6355fbbcb3a7b832ab0f62663400e6cc1471e764
--- /dev/null
+++ b/functions/basis_functions/bf_linear.h
@@ -0,0 +1,18 @@
+#ifndef BF_LINEAR_H
+#define BF_LINEAR_H
+#include "bf_base.h"
+
+/** Linear basis function */
+struct BF_Linear: public virtual BF_Base
+{
+ BF_Linear();
+ BF_Linear(const Config &c);
+ std::string get_label();
+ std::string label = "BF_Linear";
+ double epredict(const t_type &weights, const aed_rctype& aed);
+ double fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi, const size_t k);
+ force_type fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi);
+};
+#endif
diff --git a/functions/basis_functions/bf_polynomial2.cpp b/functions/basis_functions/bf_polynomial2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..be94802d700a99ea3b0e0771f216e07a2a67964b
--- /dev/null
+++ b/functions/basis_functions/bf_polynomial2.cpp
@@ -0,0 +1,52 @@
+#include "bf_polynomial2.h"
+
+Registry<Function_Base>::Register<BF_Polynomial2> BF_Polynomial2_1( "BF_Polynomial2" );
+Registry<Function_Base,Config&>::Register<BF_Polynomial2> BF_Polynomial2_2( "BF_Polynomial2" );
+
+BF_Polynomial2::BF_Polynomial2() {}
+BF_Polynomial2::BF_Polynomial2(const Config &c): Function_Base(c)
+{}
+std::string BF_Polynomial2::get_label()
+{
+ return label;
+}
+double BF_Polynomial2::epredict(const t_type &weights, const aed_rctype& aed)
+{
+ size_t b=0;
+ double res=0.0;
+ for (long int i=0; i<aed.size(); ++i) {
+ for (long int ii=i; ii<aed.size(); ++ii) {
+ res+= weights(b++)*aed(i)*aed(ii);
+ }
+ }
+ return res;
+}
+double BF_Polynomial2::fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi, const size_t k)
+{
+ double res=0.0;
+ size_t b=0;
+ for (long int i=0; i<fdij.rows(); ++i) {
+ for (long int ii=i; ii<fdij.rows(); ++ii) {
+ res -= weights(b)*(fdij(i,k)*aedi(ii) + fdij(ii,k)*aedi(i));
+ b++;
+ }
+ }
+ return res;
+}
+force_type BF_Polynomial2::fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi)
+{
+ force_type v(3);
+ v.setZero();
+ size_t b=0;
+ for (long int i=0; i<fdij.rows(); ++i) {
+ for (long int ii=i; ii<fdij.rows(); ++ii) {
+ for (size_t k=0; k<3; ++k) {
+ v(k) -= weights(b)*(fdij(i,k)*aedi(ii) + fdij(ii,k)*aedi(i));
+ }
+ b++;
+ }
+ }
+ return v;
+}
diff --git a/functions/basis_functions/bf_polynomial2.h b/functions/basis_functions/bf_polynomial2.h
new file mode 100644
index 0000000000000000000000000000000000000000..65f421292f98dd45b561dcefc26ab60091aea188
--- /dev/null
+++ b/functions/basis_functions/bf_polynomial2.h
@@ -0,0 +1,18 @@
+#ifndef BF_POLYNOMIAL2_H
+#define BF_POLYNOMIAL2_H
+#include "bf_base.h"
+
+/** Polynomial order 2 basis function with c=0 */
+struct BF_Polynomial2: public virtual BF_Base
+{
+ BF_Polynomial2();
+ BF_Polynomial2(const Config &c);
+ std::string label = "BF_Polynomial2";
+ std::string get_label();
+ double epredict(const t_type &weights, const aed_rctype& aed);
+ double fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi, const size_t k);
+ force_type fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi);
+};
+#endif
diff --git a/functions/function_base.cpp b/functions/function_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..20f43072a83025e51248f352e59ef08d56d08874
--- /dev/null
+++ b/functions/function_base.cpp
@@ -0,0 +1,2 @@
+#include "function_base.h"
+
diff --git a/functions/function_base.h b/functions/function_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..66d5823e30b56bd494aa2974d040d8f3879e1fbc
--- /dev/null
+++ b/functions/function_base.h
@@ -0,0 +1,49 @@
+#ifndef FUNCTION_BASE_h
+#define FUNCTION_BASE_h
+
+#include "../../CORE/registry.h"
+#include "../../CORE/typedefs.h"
+#include "../../CORE/config/config.h"
+
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+/** Base class for Kernels and Basis Functions */
+struct Function_Base {
+
+ private:
+ int verbose;
+
+ public:
+ Function_Base() {
+ verbose=0;
+ }
+ Function_Base(const Config &c) {
+ verbose = c.get<int>("VERBOSE");
+ }
+ void set_verbose(int v) { verbose=v; }
+ int get_verbose() { return verbose; }
+
+ // Derived classes must implement Derived() and Derived(Config)
+ virtual ~Function_Base() {};
+
+ // For Kernels
+ virtual double operator() (const aed_rctype& , const aed_rctype& ) { return 0.0; };
+ virtual aed_type derivative(const aed_rctype& , const aed_rctype& ) { return aed_type(); };
+ virtual double prime(const aed_rctype& , const aed_rctype& ,
+ const aed_rctype& ) { return 0.0; }
+ virtual void set_basis(const aeds_type ) {};
+
+ /** \brief Return label of this object */
+ virtual std::string get_label()=0;
+ virtual double epredict(const t_type &, const aed_rctype& )=0;
+ virtual double fpredict(const t_type &, const fd_type &,
+ const aed_rctype&, const size_t )=0;
+ virtual force_type fpredict(const t_type &, const fd_type &,
+ const aed_rctype& )=0;
+};
+template<> inline Registry<Function_Base>::Map Registry<Function_Base>::registry{};
+template<> inline Registry<Function_Base,Config&>::Map Registry<Function_Base,Config&>::registry{};
+#endif
diff --git a/functions/kernels/kern_all.h b/functions/kernels/kern_all.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a8704ee013a3cac753d30f69513ee792c289148
--- /dev/null
+++ b/functions/kernels/kern_all.h
@@ -0,0 +1,7 @@
+#include "kern_base.h"
+#include "kern_linear.h"
+#include "kern_rbf.h"
+#include "kern_lq.h"
+#include "kern_polynomial.h"
+#include "kern_sigmoid.h"
+#include "kern_quadratic.h"
diff --git a/functions/kernels/kern_base.cpp b/functions/kernels/kern_base.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..71791559c61fba8b3f417cd77a2a1bfaff8cd2e6
--- /dev/null
+++ b/functions/kernels/kern_base.cpp
@@ -0,0 +1,36 @@
+#include "kern_base.h"
+
+Kern_Base::~Kern_Base() {}
+void Kern_Base::set_basis(const aeds_type b)
+{
+ basis = b;
+}
+double Kern_Base::epredict(const t_type &kweights, const aed_rctype &aed)
+{
+ double energy=0;
+ for (long int b=0; b<basis.cols(); ++b) {
+ energy += kweights[b]*(*this)(basis.col(b),aed);
+ }
+ return energy;
+}
+double Kern_Base::fpredict(const t_type &kweights, const fd_type &fdij,
+ const aed_rctype &aedi, const size_t k)
+{
+ double res=0.0;
+ for (long int b=0; b<basis.cols(); ++b) {
+ res -= kweights[b]*(*this).prime(basis.col(b),aedi,fdij.col(k));
+ }
+ return res;
+}
+force_type Kern_Base::fpredict(const t_type &kweights, const fd_type &fdij,
+ const aed_rctype &aedi)
+{
+ force_type v(3);
+ v.setZero();
+ for (long int b=0; b<basis.cols(); ++b) {
+ for (size_t k=0; k<3; ++k) {
+ v(k) -= kweights[b]*(*this).prime(basis.col(b),aedi,fdij.col(k));
+ }
+ }
+ return v;
+}
diff --git a/functions/kernels/kern_base.h b/functions/kernels/kern_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..4a620487860edda9f54037f46363b1195f2f2a35
--- /dev/null
+++ b/functions/kernels/kern_base.h
@@ -0,0 +1,58 @@
+#ifndef KERN_BASE_H
+#define KERN_BASE_H
+
+#include "../../../CORE/typedefs.h"
+#include "../function_base.h"
+
+#include <iostream>
+
+/** \brief Abstract class to be used as a base for all kernels.
+ *
+ * - b = basis vector
+ * - af = atomic energy descriptor
+ * - ff = force descriptor
+ * - all derivatives are defined wrt to the second argument
+ */
+class Kern_Base: public virtual Function_Base {
+ public:
+ aeds_type basis;
+ virtual ~Kern_Base();
+
+ /**
+ * Calculate kernel value given two vectors
+ *
+ * - b = basis vector
+ * - af = atomic energy descriptor
+ */
+ virtual double operator() (const aed_rctype &b, const aed_rctype &af)=0;
+
+ /**
+ * Calculate the kernel derivative wrt to the second argument
+ *
+ * - b = basis vector
+ * - af = atomic energy descriptor
+ */
+ virtual aed_type derivative(const aed_rctype &b, const aed_rctype &af)=0;
+
+ /**
+ * Calculate inner product of the kernel derivative
+ * wrt to the second argument (b wrt af) with the force descriptor
+ *
+ * - b = basis vector
+ * - af = atomic energy descriptor
+ * - ff = force descriptor (TODO i-th dir component of it)
+ * TODO consider calculating all 3-directions at once
+ */
+ virtual double prime(const aed_rctype &b, const aed_rctype &af,
+ const aed_rctype &ff)=0;
+
+ /** \brief Set basis for calculations */
+ virtual void set_basis(const aeds_type b);
+
+ virtual double epredict(const t_type &kweights, const aed_rctype &aed);
+ virtual double fpredict(const t_type &kweights, const fd_type &fdij,
+ const aed_rctype &aedi, const size_t k);
+ virtual force_type fpredict(const t_type &kweights, const fd_type &fdij,
+ const aed_rctype &aedi);
+};
+#endif
diff --git a/functions/kernels/kern_linear.cpp b/functions/kernels/kern_linear.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3c80fb7579e4a0b88aae38e5e3ab8d7c8b41c5d0
--- /dev/null
+++ b/functions/kernels/kern_linear.cpp
@@ -0,0 +1,40 @@
+#include "kern_linear.h"
+
+Registry<Function_Base>::Register<Kern_Linear> Kern_Linear_1( "Kern_Linear" );
+Registry<Function_Base,Config&>::Register<Kern_Linear> Kern_Linear_2( "Kern_Linear" );
+
+Kern_Linear::Kern_Linear() {}
+Kern_Linear::Kern_Linear (const Config &c):
+ Function_Base(c)
+{}
+double Kern_Linear::operator() (const aed_rctype& b, const aed_rctype& af)
+{
+ return b.transpose()*af;;
+}
+aed_type Kern_Linear::derivative(const aed_rctype& b, const aed_rctype&)
+{
+ return b;
+}
+double Kern_Linear::prime(const aed_rctype& b, const aed_rctype& , const aed_rctype& ff)
+{
+ // return derivative(b,af).transpose() * ff;
+ return b.transpose() * ff;
+}
+std::string Kern_Linear::get_label()
+{
+ return label;
+}
+double Kern_Linear::epredict(const t_type &weights, const aed_rctype& aed)
+{
+ return weights.transpose()*aed;
+}
+double Kern_Linear::fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& , const size_t k)
+{
+ return -fdij.col(k).transpose() * weights;
+}
+force_type Kern_Linear::fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& )
+{
+ return -fdij.transpose() * weights;
+}
diff --git a/functions/kernels/kern_linear.h b/functions/kernels/kern_linear.h
new file mode 100644
index 0000000000000000000000000000000000000000..cbad9db39c5f008bd10ff61b5bfe7e966143488d
--- /dev/null
+++ b/functions/kernels/kern_linear.h
@@ -0,0 +1,36 @@
+#ifndef KERN_LINEAR_H
+#define KERN_LINEAR_H
+#include "kern_base.h"
+/**
+ * Linear kernel also knows as dot product kernel
+ *
+ * Defined for two vectors **x** and **y**:
+ *
+ * \f[
+ * K(\mathbf{x}, \mathbf{y}) = \mathbf{x}^T \mathbf{y} = \sum_i x^{(i)} y^{(i)}
+ * \f]
+ *
+ * @see Kern_Base BF_Linear
+ */
+class Kern_Linear : public virtual Kern_Base {
+ public:
+ Kern_Linear ();
+ Kern_Linear (const Config &c);
+ std::string get_label();
+
+ /**
+ * Label used for this class
+ */
+ std::string label = "Kern_Linear";
+ double operator() (const aed_rctype& b, const aed_rctype& af);
+ aed_type derivative(const aed_rctype& b, const aed_rctype& );
+ double prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff);
+ void set_basis(const aeds_type ) {}
+
+ double epredict(const t_type &weights, const aed_rctype& aed);
+ double fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi, const size_t k);
+ force_type fpredict(const t_type &weights, const fd_type &fdij,
+ const aed_rctype& aedi);
+};
+#endif
diff --git a/functions/kernels/kern_lq.cpp b/functions/kernels/kern_lq.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..82b3db722d0667921886de8f3dc1bf88f2e66ac3
--- /dev/null
+++ b/functions/kernels/kern_lq.cpp
@@ -0,0 +1,29 @@
+#include "kern_lq.h"
+
+Registry<Function_Base>::Register<Kern_LQ> Kern_LQ_1( "Kern_LQ" );
+Registry<Function_Base,Config&>::Register<Kern_LQ> Kern_LQ_2( "Kern_LQ" );
+
+Kern_LQ::Kern_LQ () {}
+Kern_LQ::Kern_LQ (const Config &c):
+ Function_Base(c)
+{
+ if (get_verbose()) std::cout << std::endl << label << std::endl;
+}
+double Kern_LQ::operator() (const aed_rctype& b, const aed_rctype& af)
+{
+ double x = b.transpose()*af;
+ return x*x + x;
+}
+aed_type Kern_LQ::derivative(const aed_rctype& b, const aed_rctype& af)
+{
+ double temp = 2.0 * b.transpose()*af;
+ return temp*b + b;
+}
+double Kern_LQ::prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff)
+{
+ return derivative(b, af).transpose()*ff;
+}
+std::string Kern_LQ::get_label()
+{
+ return label;
+}
diff --git a/functions/kernels/kern_lq.h b/functions/kernels/kern_lq.h
new file mode 100644
index 0000000000000000000000000000000000000000..9d5bd19b3ffb3684715b73c9865f3b03f042bffe
--- /dev/null
+++ b/functions/kernels/kern_lq.h
@@ -0,0 +1,32 @@
+#ifndef KERN_LQ_H
+#define KERN_LQ_H
+
+#include "kern_base.h"
+/**
+ * Linear + Quadratic kernel
+ *
+ * Defined for two vectors **x** and **y**:
+ *
+ * \f[
+ * K(\mathbf{x}, \mathbf{y}) = \mathbf{x}^T \mathbf{y} +
+ * \Big(\mathbf{x}^T \mathbf{y} \Big)^2
+ * = \sum_i x^{(i)} y^{(i)} + \Big(\sum_i x^{(i)} y^{(i)} \Big)^2
+ * \f]
+ *
+ * @see Kern_Base Kern_Linear Kern_Quadratic
+ */
+class Kern_LQ : public virtual Kern_Base {
+ public:
+ Kern_LQ ();
+ Kern_LQ (const Config &c);
+ /**
+ * Label used for this class
+ */
+ std::string label = "Kern_LQ";
+ double operator() (const aed_rctype& b, const aed_rctype& af);
+ aed_type derivative(const aed_rctype& b, const aed_rctype& af);
+ double prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff);
+ std::string get_label();
+
+};
+#endif
diff --git a/functions/kernels/kern_polynomial.cpp b/functions/kernels/kern_polynomial.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..01dcf5e21f8f69e3ae29f414211881b405c1beb7
--- /dev/null
+++ b/functions/kernels/kern_polynomial.cpp
@@ -0,0 +1,31 @@
+#include "kern_polynomial.h"
+
+Registry<Function_Base>::Register<Kern_Polynomial> Kern_Polynomial_1( "Kern_Polynomial" );
+Registry<Function_Base,Config&>::Register<Kern_Polynomial> Kern_Polynomial_2( "Kern_Polynomial" );
+
+Kern_Polynomial::Kern_Polynomial () {}
+Kern_Polynomial::Kern_Polynomial (const Config &c):
+ Function_Base(c),
+ d(c.get<int>("MPARAMS",0)),
+ gamma(c.get<double>("MPARAMS",1)),
+ c(c.get<double>("MPARAMS",2))
+{
+ if (get_verbose()) std::cout << std::endl << label << " | degree: " << d
+ <<" | gamma: " << gamma << " | c: " << c << std::endl;
+}
+double Kern_Polynomial::operator() (const aed_rctype& b, const aed_rctype& af)
+{
+ return std::pow(gamma*(b.dot(af))+c,d);
+}
+aed_type Kern_Polynomial::derivative(const aed_rctype& b, const aed_rctype& af)
+{
+ return d*gamma*std::pow(gamma*(b.dot(af))+c,d-1)*b;
+}
+double Kern_Polynomial::prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff)
+{
+ return derivative(b, af).transpose()*ff;
+}
+std::string Kern_Polynomial::get_label()
+{
+ return label;
+}
diff --git a/functions/kernels/kern_polynomial.h b/functions/kernels/kern_polynomial.h
new file mode 100644
index 0000000000000000000000000000000000000000..02bf26262d3d044c338a2e1fc7e1e9a8b4d853b9
--- /dev/null
+++ b/functions/kernels/kern_polynomial.h
@@ -0,0 +1,34 @@
+#ifndef KERN_POLYNOMIAL_H
+#define KERN_POLYNOMIAL_H
+#include "kern_base.h"
+/** \brief Polynomial kernel.
+ *
+ * Defined for two vectors **x** and **y**:
+ *
+ * \f[
+ * K(\mathbf{x}, \mathbf{y}) = \big( \gamma*\mathbf{x}^T \mathbf{y} +
+ * \mathrm{c} \big)^{\mathrm{d}}
+ * \f]
+ *
+ * Required Config key: \ref MPARAMS <int> d <double> \f$ \gamma \f$ c
+ *
+ * @see Kern_Base
+ */
+class Kern_Polynomial : public virtual Kern_Base {
+ public:
+ Kern_Polynomial ();
+ Kern_Polynomial (const Config &c);
+ int d; // power
+ double gamma; // scalling
+ double c; // shift
+ /**
+ * Label used for this class
+ */
+ std::string label = "Kern_Polynomial";
+ std::string get_label();
+ double operator() (const aed_rctype& b, const aed_rctype& af);
+ aed_type derivative(const aed_rctype& b, const aed_rctype& af);
+ double prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff);
+
+};
+#endif
diff --git a/functions/kernels/kern_quadratic.cpp b/functions/kernels/kern_quadratic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0d61f734710c1a39e86fa23ab88d2c9a5858a214
--- /dev/null
+++ b/functions/kernels/kern_quadratic.cpp
@@ -0,0 +1,32 @@
+#include "kern_quadratic.h"
+
+Registry<Function_Base>::Register<Kern_Quadratic> Kern_Quadratic_1( "Kern_Quadratic" );
+Registry<Function_Base,Config&>::Register<Kern_Quadratic> Kern_Quadratic_2( "Kern_Quadratic" );
+
+Kern_Quadratic::Kern_Quadratic () {}
+Kern_Quadratic::Kern_Quadratic (const Config &c):
+ Function_Base(c)
+{
+ if (get_verbose()) std::cout << std::endl << label << std::endl;
+}
+double Kern_Quadratic::operator() (const aed_rctype& b, const aed_rctype& af)
+{
+ double x = b.transpose()*af;
+ return x*x;
+}
+aed_type Kern_Quadratic::derivative(const aed_rctype& b, const aed_rctype& af)
+{
+ double temp = 2.0 * b.transpose()*af;
+ return b*temp;;
+}
+double Kern_Quadratic::prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff)
+{
+ //return dot(derivative(b, af), ff);
+ // it is faster this way
+ double scale = 2.0*b.transpose()*af;
+ return scale*ff.transpose()*b;
+}
+std::string Kern_Quadratic::get_label()
+{
+ return label;
+}
diff --git a/functions/kernels/kern_quadratic.h b/functions/kernels/kern_quadratic.h
new file mode 100644
index 0000000000000000000000000000000000000000..d9e886bb54253afc4ee918bf98319f117887cc45
--- /dev/null
+++ b/functions/kernels/kern_quadratic.h
@@ -0,0 +1,30 @@
+#ifndef KERN_QUADRATIC_H
+#define KERN_QUADRATIC_H
+#include "kern_base.h"
+/**
+ * Quadratic kernel - special case of 2nd order polynomial kernel
+ *
+ * Defined for two vectors **x** and **y**:
+ *
+ * \f[
+ * K(\mathbf{x}, \mathbf{y}) = \Big( \mathbf{x}^T \mathbf{y} \Big)^2
+ * = \Big( \sum_i x^{(i)} y^{(i)} \Big)^2
+ * \f]
+ *
+ * @see Kern_Base
+ */
+class Kern_Quadratic : public virtual Kern_Base {
+ public:
+ Kern_Quadratic ();
+ Kern_Quadratic (const Config &c);
+ /**
+ * Label used for this class
+ */
+ std::string label = "Kern_Quadratic";
+ double operator() (const aed_rctype& b, const aed_rctype& af);
+ aed_type derivative(const aed_rctype& b, const aed_rctype& af);
+ double prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff);
+ std::string get_label();
+
+};
+#endif
diff --git a/functions/kernels/kern_rbf.cpp b/functions/kernels/kern_rbf.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a96a9352feb4e8462f1f5ad8fd1eb6c4428ad93d
--- /dev/null
+++ b/functions/kernels/kern_rbf.cpp
@@ -0,0 +1,31 @@
+#include "kern_rbf.h"
+
+Registry<Function_Base>::Register<Kern_RBF> Kern_RBF_1( "Kern_RBF" );
+Registry<Function_Base,Config&>::Register<Kern_RBF> Kern_RBF_2( "Kern_RBF" );
+
+Kern_RBF::Kern_RBF () {}
+Kern_RBF::Kern_RBF (const Config &c):
+ Function_Base(c),
+ sigma(c.get<double>("MPARAMS")),
+ gamma(1.0/(2.0*sigma*sigma))
+{
+ if (get_verbose()) std::cout << std::endl << label << " | sigma: " << sigma
+ << " | gamma: " << gamma << std::endl;
+}
+double Kern_RBF::operator() (const aed_rctype& b, const aed_rctype& af)
+{
+ const aed_type v = af-b;
+ return exp(-gamma*(v.dot(v)));
+}
+aed_type Kern_RBF::derivative(const aed_rctype& b, const aed_rctype& af)
+{
+ return -2.0*gamma*(af-b)*(*this)(af,b);
+}
+double Kern_RBF::prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff)
+{
+ return derivative(b, af).transpose()*ff;
+}
+std::string Kern_RBF::get_label()
+{
+ return label;
+}
diff --git a/functions/kernels/kern_rbf.h b/functions/kernels/kern_rbf.h
new file mode 100644
index 0000000000000000000000000000000000000000..d91c493ffa0ed3a5f89bc17471968af41dc40419
--- /dev/null
+++ b/functions/kernels/kern_rbf.h
@@ -0,0 +1,35 @@
+#ifndef KERN_RBF_H
+#define KERN_RBF_H
+#include "kern_base.h"
+/** \brief Radial Basis Function kernel
+ *
+ * Defined for two vectors **x** and **y**:
+ *
+ * \f[
+ * K(\mathbf{x}, \mathbf{y}) = \exp\Big( -\frac{||\mathbf{x}-\mathbf{y}||^2}
+ * {2\sigma^2} \Big)
+ * \f]
+ *
+ * Required Config key: \ref MPARAMS <double> \f$ \sigma \f$
+ *
+ * @see Kern_Base Kern_Linear Kern_Quadratic
+ */
+class Kern_RBF : public virtual Kern_Base {
+ public:
+ double sigma;
+ /**
+ * Label used for this class
+ */
+ Kern_RBF ();
+ Kern_RBF (const Config &c);
+ std::string label = "Kern_RBF";
+ double operator() (const aed_rctype& b, const aed_rctype& af);
+ aed_type derivative(const aed_rctype& b, const aed_rctype& af);
+ double prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff);
+ std::string get_label();
+
+ private:
+ double gamma;
+
+};
+#endif
diff --git a/functions/kernels/kern_sigmoid.cpp b/functions/kernels/kern_sigmoid.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6f0c9470844d9b359ff2394250836db03d19f1f2
--- /dev/null
+++ b/functions/kernels/kern_sigmoid.cpp
@@ -0,0 +1,30 @@
+#include "kern_sigmoid.h"
+
+Registry<Function_Base>::Register<Kern_Sigmoid> Kern_Sigmoid_1( "Kern_Sigmoid" );
+Registry<Function_Base,Config&>::Register<Kern_Sigmoid> Kern_Sigmoid_2( "Kern_Sigmoid" );
+
+Kern_Sigmoid::Kern_Sigmoid () {}
+Kern_Sigmoid::Kern_Sigmoid (const Config &c):
+ Function_Base(c),
+ gamma(c.get<double>("MPARAMS",0)),
+ c(c.get<double>("MPARAMS",1))
+{
+ if (get_verbose()) std::cout << std::endl << label << " | gamma: " << gamma
+ << " | c: " << c << std::endl;
+}
+double Kern_Sigmoid::operator() (const aed_rctype& b, const aed_rctype& af)
+{
+ return std::tanh(gamma*(b.dot(af))+c);
+}
+aed_type Kern_Sigmoid::derivative(const aed_rctype& b, const aed_rctype& af)
+{
+ return gamma*(1.0-std::pow((*this)(af,b),2))*b;
+}
+double Kern_Sigmoid::prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff)
+{
+ return derivative(b, af).transpose()*ff;
+}
+std::string Kern_Sigmoid::get_label()
+{
+ return label;
+}
diff --git a/functions/kernels/kern_sigmoid.h b/functions/kernels/kern_sigmoid.h
new file mode 100644
index 0000000000000000000000000000000000000000..600251af6dee361421fc755d841b72ad45033165
--- /dev/null
+++ b/functions/kernels/kern_sigmoid.h
@@ -0,0 +1,33 @@
+#ifndef KERN_SIGMOID_H
+#define KERN_SIGMOID_H
+#include "kern_base.h"
+/** \brief Sigmoid kernel.
+ *
+ * Defined for two vectors **x** and **y**:
+ *
+ * \f[
+ * K(\mathbf{x}, \mathbf{y}) = \tanh\big( \gamma*\mathbf{x}^T \mathbf{y}
+ * + \mathrm{c} \big)
+ * \f]
+ *
+ * Required Config key: \ref MPARAMS <double> \f$ \gamma \f$ c
+ *
+ * @see Kern_Base
+ */
+class Kern_Sigmoid : public virtual Kern_Base {
+ public:
+ double gamma; // scalling
+ double c; // shift
+ /**
+ * Label used for this class
+ */
+ Kern_Sigmoid ();
+ Kern_Sigmoid (const Config &c);
+ std::string label = "Kern_Sigmoid";
+ double operator() (const aed_rctype& b, const aed_rctype& af);
+ aed_type derivative(const aed_rctype& b, const aed_rctype& af);
+ double prime(const aed_rctype& b, const aed_rctype& af, const aed_rctype& ff);
+ std::string get_label();
+
+};
+#endif
diff --git a/linear_regressor.h b/linear_regressor.h
new file mode 100644
index 0000000000000000000000000000000000000000..feb69fdd5f6d9b984f3cb5b9402c9e68fc9e1b34
--- /dev/null
+++ b/linear_regressor.h
@@ -0,0 +1,111 @@
+#ifndef LINEAR_REGRESSOR_H
+#define LINEAR_REGRESSOR_H
+
+#include "../CORE/config/config.h"
+#include "../CORE/typedefs.h"
+#include "../LIBS/Eigen/Dense"
+#include "ea.h"
+#include <vector>
+class LinearRegressor {
+ // do training...
+ // In general there are 3 methods to solve least squares problem:
+ // 1) Normal equation -direct inverse - fast but least accurate:
+ //Eigen::MatrixXd I = Eigen::MatrixXd::Identity(Phi.cols(),Phi.cols());
+ //Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> PhiPhi = Phi.transpose()*Phi;
+ //weights = (PhiPhi+lambda*I).inverse() * Phi.transpose() * T;
+ // 2) QR - in between 1) and 3) in terms of speed and accuracy
+ //weights = Phi.colPivHouseholderQr().solve(T);
+ // 3) SVD - most accurate but the slowest
+ // w = V*D*S(*)*U^T*T
+ // where S(*) = S.asDiagonal().inverse()
+ // and D is a filter matrix with diagonal entires: S(i)^2/(S(i)^2+lambda)
+ public:
+ static void train(Config &config, const phi_type &Phi, const t_type &T,
+ t_type &weights, mat &Sigma) {
+ int verbose = config.get<int>("VERBOSE");
+
+ double lambda=config.get<double>("LAMBDA");
+
+ if (lambda == 0)
+ weights = Phi.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(T);
+ else {
+ double alpha = config.get<double>("ALPHA");
+ double beta = config.get<double>("BETA");
+ //Eigen::JacobiSVD<mat> svd(Phi, Eigen::ComputeThinU | Eigen::ComputeThinV);
+ Eigen::BDCSVD<mat> svd(Phi, Eigen::ComputeThinU | Eigen::ComputeThinV);
+ mat S =svd.singularValues();
+ const mat &U =svd.matrixU();
+ const mat &V =svd.matrixV();
+
+ // Filter out v.small singularvalues from S and its inverse
+ for (long int i=0; i<S.rows(); ++i) {
+ if (S(i,0)<std::numeric_limits<double>::min()) {
+ S(i,0) = 0.0;
+ }
+ }
+
+ if (lambda<0) {
+ EA ea(config);
+ ea.run(Phi,T,S,U,V,alpha,beta);
+ lambda=alpha/beta;
+
+ mat P_inv(S.rows(),S.rows());
+ P_inv.setZero();
+
+ for (long int i=0; i<S.rows(); ++i) {
+ double s = S(i,0);
+ P_inv(i,i) = 1.0/(beta*s*s+alpha);
+ }
+
+ Sigma = V*P_inv*V.transpose();
+ config.remove("ALPHA");
+ config.remove("BETA");
+ config.add("ALPHA",alpha);
+ config.add("BETA",beta);
+ }
+ mat D(S.rows(),S.rows());
+ D.setZero();
+ for (long int i=0; i<S.rows(); ++i) {
+ double s = S(i,0);
+ if (s>std::numeric_limits<double>::min())
+ D(i,i) = (s)/(s*s+lambda);
+ else
+ D(i,i) = 0.0;
+ }
+ if (verbose) std::cout << std::endl << "REG LAMBDA: " << lambda << std::endl;
+ weights = V*D*U.transpose()*T;
+ //weights = beta*Sigma*V*S.asDiagonal()*U.transpose()*T;
+ // add Sigma to the config file
+ dump_sigma(config,Sigma);
+ }
+ }
+ static void dump_sigma(Config &config, mat &Sigma) {
+ vec Sigma_as_vector(Eigen::Map<vec>(Sigma.data(), Sigma.cols()*Sigma.rows()));
+ config.add("SIGMA",Sigma.rows());
+ config.add<vec>("SIGMA",Sigma_as_vector);
+
+ }
+ static void read_sigma(Config &config, mat &Sigma) {
+ using vec=std::vector<double>;
+ // get N = number of rows or cols, Sigma is NxN
+ size_t N;
+ try {
+ N = config.get<size_t>("SIGMA");
+ }
+ catch (const std::runtime_error& error) {
+ throw std::runtime_error("Sigma matrix is not computed. \n\
+ Hint: It is only computed for LAMBDA -1.\n");
+ }
+ // +1 as first number is N - a number of rows and cols
+ vec Sigma_as_vector(N*N+1);
+ // read sigma including N, it is sotred col by col
+ config.get<vec>("SIGMA",Sigma_as_vector);
+ Sigma_as_vector.erase(Sigma_as_vector.begin());
+ Sigma.resize(N,N);
+ for (size_t c=0; c<N; ++c)
+ for (size_t r=0; r<N; ++r)
+ Sigma(r,c) = Sigma_as_vector.at(N*c+r);
+
+ }
+};
+#endif
diff --git a/m_base.h b/m_base.h
new file mode 100644
index 0000000000000000000000000000000000000000..4ecb44b1d7a415c2795d466a0f61fec80c5a6248
--- /dev/null
+++ b/m_base.h
@@ -0,0 +1,35 @@
+#ifndef M_BASE_H
+#define M_BASE_H
+
+#include "../LIBS/Eigen/Dense"
+#include "../CORE/typedefs.h"
+#include "functions/function_base.h"
+
+class M_Base {
+
+ public:
+ //Normaliser norm; // TODO?
+ virtual ~M_Base() {};
+
+ /** \brief Predict local energy of an atom or bond energy.
+ *
+ * The result depends on how aed is computed.
+ *
+ * If it is computed between a pair of atoms than the result is a bond energy.
+ *
+ * If aed contains sum over all nearest neighbours than the result is
+ * a local atomic energy \f$ E_i \f$.
+ * */
+ virtual double epredict(const aed_rctype &aed)=0;
+
+ /** \brief Predict force between a pair of atoms in a k-direction. */
+ virtual double fpredict(const fd_type &fdij, const aed_rctype &aedi, size_t k)=0;
+
+ /** \brief Predict force between a pair of atoms. */
+ virtual force_type fpredict(const fd_type &fdij,
+ const aed_rctype &aedi)=0;
+
+
+};
+//template<> inline Registry<M_Base,Function_Base&,Config&>::Map Registry<M_Base,Function_Base&,Config&>::registry{};
+#endif
diff --git a/m_blr_core.h b/m_blr_core.h
new file mode 100644
index 0000000000000000000000000000000000000000..dea3ea2250f6ba1829a73780aca314998417a73b
--- /dev/null
+++ b/m_blr_core.h
@@ -0,0 +1,73 @@
+#ifndef M_BLR_Core_H
+#define M_BLR_Core_H
+
+#include "linear_regressor.h"
+#include "../CORE/config/config.h"
+#include "functions/function_base.h"
+#include "functions/basis_functions/bf_base.h"
+#include "m_core.h"
+
+#include <iostream>
+
+template
+<class BF=Function_Base&>
+class M_BLR_Core: public virtual M_Core {
+ public:
+ Config &config;
+ mat Sigma;
+ BF bf;
+ M_BLR_Core(Config &c):
+ config(c),
+ bf(c)
+ {
+ static_assert(std::is_base_of<BF_Base, BF>::value,
+ "\nThe provided Basis function is not of a BasisFunction type.\n\
+ For BLR use basis functions defined in basis_functions.h.\n\
+ If you want to use this kernel functions use M_KRR class instead.\n\
+ e.g. Use BF_Linear instead of Kern_Linear\n.");
+
+ init();
+ }
+
+ M_BLR_Core(BF &bf, Config &c):
+ config(c),
+ bf(bf)
+ {
+ if (dynamic_cast<BF_Base*>(&bf) == nullptr)
+ throw std::invalid_argument("Provided object is not of BF_Base type");
+
+ //static_assert(std::is_same<Function_Base&, BF>::value,
+ // "This constructor requires BF=Function_Base&\n");
+
+ init();
+ }
+ double predict(const rvec &v) const{
+ return bf.epredict(get_weights(),v);
+ };
+
+ void train(phi_type &Phi, const t_type &T) {
+ if(trained) {
+ throw std::runtime_error("This object is already trained!");
+ }
+ LinearRegressor::train(config,Phi, T,weights,Sigma);
+ trained=true;
+ }
+ t_type get_weights_uncertainty() const{
+ double lambda=config.get<double>("LAMBDA");
+ if(lambda >= 0) throw std::runtime_error(
+ "Sigma matrix is only computed for LAMBDA < 0");
+ return Sigma.diagonal();
+ }
+
+ private:
+ void init() {
+ if (config.exist("WEIGHTS")) {
+ weights.resize(config.size("WEIGHTS"));
+ config.get<t_type>("WEIGHTS",weights);
+ trained=true;
+ }
+
+ M_Core::verbose=(config.get<int>("VERBOSE"));
+ }
+};
+#endif
diff --git a/m_core.h b/m_core.h
new file mode 100644
index 0000000000000000000000000000000000000000..1caa4b00da4379be1ec5c896352c98184202d2a0
--- /dev/null
+++ b/m_core.h
@@ -0,0 +1,28 @@
+#ifndef M_Core_H
+#define M_Core_H
+
+#include "../CORE/typedefs.h"
+
+class M_Core {
+ public:
+ int verbose;
+ bool trained=false;
+ t_type weights;
+
+ virtual ~M_Core() {}
+
+ virtual double predict(const rvec &v)const=0;
+ virtual void train(phi_type &Phi, const t_type &T)=0;
+ virtual t_type get_weights_uncertainty()const=0;
+
+ bool is_trained() const {
+ return trained;
+ }
+ const t_type &get_weights() const {
+ return weights;
+ }
+ void set_weights(const t_type w) {
+ weights=w;
+ }
+};
+#endif
diff --git a/m_krr_core.h b/m_krr_core.h
new file mode 100644
index 0000000000000000000000000000000000000000..3499aad027b3ccd01c04a062c82f24e3d69f5ea6
--- /dev/null
+++ b/m_krr_core.h
@@ -0,0 +1,81 @@
+#ifndef M_KRR_Core_H
+#define M_KRR_Core_H
+
+#include "linear_regressor.h"
+#include "../CORE/config/config.h"
+#include "functions/function_base.h"
+#include "functions/kernels/kern_base.h"
+#include "m_core.h"
+#include "ekm.h"
+
+#include <iostream>
+
+template
+<class K=Function_Base&>
+class M_KRR_Core: public virtual M_Core {
+ public:
+ Config &config;
+ mat Sigma;
+ K kernel;
+ EKM<K> ekm;
+ M_KRR_Core(Config &c):
+ config(c),
+ kernel(c),
+ ekm(c)
+ {
+ static_assert(std::is_base_of<Kern_Base, K>::value,
+ "\nThe provided Kernel K is not of a Kernel type.\n\
+ For KRR use kernels defined in kernel.h.\n\
+ e.g. Use Kern_Linear instead of BF_Linear\n.");
+
+ init();
+ }
+
+ M_KRR_Core(K &kernel, Config &c):
+ config(c),
+ kernel(kernel),
+ ekm(kernel)
+ {
+ if (dynamic_cast<Kern_Base*>(&kernel) == nullptr)
+ throw std::invalid_argument("Provided object is not of Kern_Base type");
+
+ init();
+ }
+ double predict(const rvec &v) const{
+ return kernel.epredict(get_weights(),v);
+ };
+
+ void train(phi_type &Phi, const t_type &T) {
+ if(trained) {
+ throw std::runtime_error("This object is already trained!");
+ }
+ if (kernel.get_label()!="Kern_Linear") {
+ ekm.project(Phi);
+ }
+ LinearRegressor::train(config,Phi, T,weights,Sigma);
+
+ if (kernel.get_label()!="Kern_Linear") {
+ //kernalize weights
+ weights = ekm.KK.transpose()*weights;
+ }
+ trained=true;
+ }
+ t_type get_weights_uncertainty() const{
+ double lambda=config.get<double>("LAMBDA");
+ if(lambda >= 0) throw std::runtime_error(
+ "Sigma matrix is only computed for LAMBDA < 0");
+ return Sigma.diagonal();
+ }
+
+ private:
+ void init() {
+ if (config.exist("WEIGHTS")) {
+ weights.resize(config.size("WEIGHTS"));
+ config.get<t_type>("WEIGHTS",weights);
+ trained=true;
+ }
+
+ M_Core::verbose=(config.get<int>("VERBOSE"));
+ }
+};
+#endif