Init

cutoffs/cut_all.h

+#include "cutoffs.h"

cutoffs/cutoffs.cpp

+#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;
+}

cutoffs/cutoffs.h

+#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

dc_selector.h

+#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

descriptors/d2/d2_all.h

+#include "d2_base.h"
+#include "d2_dummy.h"
+#include "d2_blip.h"
+#include "d2_bp.h"
+#include "d2_eam.h"
+#include "d2_lj.h"

descriptors/d2/d2_base.cpp

+#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);
+}

descriptors/d2/d2_base.h

+#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

descriptors/d2/d2_blip.cpp

+#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;
+}

descriptors/d2/d2_blip.h

+#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

descriptors/d2/d2_bp.cpp

+#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;
+}

descriptors/d2/d2_bp.h

+#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

descriptors/d2/d2_dummy.cpp

+#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; }

descriptors/d2/d2_dummy.h

+#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

descriptors/d2/d2_eam.cpp

+#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;
+	}
+}
+

descriptors/d2/d2_eam.h

+#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

descriptors/d2/d2_lj.cpp

+#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;
+}

descriptors/d2/d2_lj.h

+#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

descriptors/d2/d2_mie.cpp

+#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;
+}

descriptors/d2/d2_mie.h

+#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

descriptors/d3/d3_all.h

+#include "d3_base.h"
+#include "d3_dummy.h"