diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0558d2011455be7bb71d95f66210afa4118481ae..7e5a64d392231dcc7d5b0de3b162f55230704445 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -68,7 +68,7 @@ endif()
message("${TADAH} Develop build: Directory not found.")
endif()
endfunction()
- find_directory_case_insensitive("models" "$ENV{TADAH_PATH}" DEVELOP_CORE_DIR)
+ find_directory_case_insensitive("models" "$ENV{TADAH_PATH}" DEVELOP_MODELS_DIR)
set(FETCHCONTENT_SOURCE_DIR_TADAH.MODELS "${DEVELOP_MODELS_DIR}")
endif()
diff --git a/include/tadah/mlip/analytics/analytics.h b/include/tadah/mlip/analytics/analytics.h
index 71cda1da7ecc506fa2e8869e8cfe41ce2c34a1ff..5c172f9435b98fd9acb9bb0292409eff235b90c5 100644
--- a/include/tadah/mlip/analytics/analytics.h
+++ b/include/tadah/mlip/analytics/analytics.h
@@ -19,40 +19,40 @@ class Analytics {
Analytics(const StructureDB &st, const StructureDB &stp);
/** Return Energy/atom Mean Absolute Error for each DBFILE. */
- t_type calc_e_mae() const;
+ tadah::core::t_type calc_e_mae() const;
/** Return Force Mean Absolute Error for each DBFILE. */
- t_type calc_f_mae() const;
+ tadah::core::t_type calc_f_mae() const;
/** Return Stress Mean Absolute Error for each DBFILE.
*
* Calculated using 6 independed components.
*/
- t_type calc_s_mae() const;
+ tadah::core::t_type calc_s_mae() const;
/** Return Energy/atom Relative Root Mean Square Error for each DBFILE. */
- t_type calc_e_rrmse() const;
+ tadah::core::t_type calc_e_rrmse() const;
/** Return Energy/atom Root Mean Square Error for each DBFILE. */
- t_type calc_e_rmse() const;
+ tadah::core::t_type calc_e_rmse() const;
/** Return Force Root Mean Square Error for each DBFILE. */
- t_type calc_f_rmse() const;
+ tadah::core::t_type calc_f_rmse() const;
/** Return Stress Root Mean Square Error for each DBFILE.
*
* Calculated using 6 independed components.
*/
- t_type calc_s_rmse() const;
+ tadah::core::t_type calc_s_rmse() const;
/** Return Energy/atom coefficient of determination for each DBFILE. */
- t_type calc_e_r_sq() const;
+ tadah::core::t_type calc_e_r_sq() const;
/** Return Force coefficient of determination for each DBFILE. */
- t_type calc_f_r_sq() const;
+ tadah::core::t_type calc_f_r_sq() const;
/** Return Stress coefficient of determination for each DBFILE. */
- t_type calc_s_r_sq() const;
+ tadah::core::t_type calc_s_r_sq() const;
};
}
diff --git a/include/tadah/mlip/analytics/statistics.h b/include/tadah/mlip/analytics/statistics.h
index 4bd75bfc92cbccaf9b224c0666da0a62228544f4..c48bea09a5421c49dd204f5ba498b1e26edb3548 100644
--- a/include/tadah/mlip/analytics/statistics.h
+++ b/include/tadah/mlip/analytics/statistics.h
@@ -3,9 +3,12 @@
#include <tadah/core/core_types.h>
+namespace tadah {
+namespace mlip {
+
/** Some basis statistical tools */
class Statistics {
- using vec = aed_type;
+ using vec = tadah::core::aed_type;
public:
/** Residual sum of squares. */
static double res_sum_sq(const vec &obs, const vec &pred);
@@ -23,4 +26,6 @@ class Statistics {
static double mean(const vec &v);
};
+}
+}
#endif
diff --git a/include/tadah/mlip/descriptors_calc.h b/include/tadah/mlip/descriptors_calc.h
index 3bb374319c8cec040f2241c64de9f5a5de911381..2e95c37efd4e0a9ff27692f0b40e2f706ab5431a 100644
--- a/include/tadah/mlip/descriptors_calc.h
+++ b/include/tadah/mlip/descriptors_calc.h
@@ -1,7 +1,6 @@
#ifndef DESCRIPTORS_CALC_H
#define DESCRIPTORS_CALC_H
-#include "tadah/mlip/structure_neighbor_view.h"
#include <tadah/mlip/descriptors_calc_base.h>
#include <tadah/mlip/structure.h>
#include <tadah/mlip/structure_db.h>
@@ -13,6 +12,8 @@
#include <tadah/models/descriptors/d3/d3_base.h>
#include <tadah/models/descriptors/dm/dm_base.h>
#include <tadah/models/cutoffs.h>
+#include <tadah/mlip/neighbor_list_db.h>
+#include <tadah/mlip/structure_neighbor_view.h>
namespace tadah {
namespace mlip {
@@ -22,7 +23,7 @@ namespace mlip {
*
* Usage example:
* \code{.cpp}
- * # Build calculator using Config object.
+ * # Build calculator using tadah::core::Config object.
* DescriptorsCalc<D2_LJ, D3_Dummy, DM_Dummy> dc(config);
*
* # Calculate all descriptors in a Structure st
@@ -33,7 +34,7 @@ namespace mlip {
* \endcode
*
* This object sets the following \ref INTERNAL_KEYS
- * to the provided Config object:
+ * to the provided tadah::core::Config object:
*
* \ref SIZE2B,\ref SIZE3B,\ref SIZEMB,\ref DSIZE
*
@@ -42,18 +43,18 @@ namespace mlip {
* \ref RCUT2B,\ref RCUT3B,\ref RCUTMB (for initialised ones)
*
*
- * @tparam D2 D2_Base child, two-body type calculator.
- * @tparam D3 D3_Base child, three-body type calculator.
- * @tparam DM DM_Base child, many-body type calculator.
- * @tparam C2 Cut_Base child, two-body cutoff function.
- * @tparam C3 Cut_Base child, three-body cutoff function.
- * @tparam CM Cut_Base child, many-body cutoff function.
+ * @tparam D2 tadah::models::D2_Base child, two-body type calculator.
+ * @tparam D3 tadah::models::D3_Base child, three-body type calculator.
+ * @tparam DM tadah::models::DM_Base child, many-body type calculator.
+ * @tparam C2 tadah::models::Cut_Base child, two-body cutoff function.
+ * @tparam C3 tadah::models::Cut_Base child, three-body cutoff function.
+ * @tparam CM tadah::models::Cut_Base child, many-body cutoff function.
*/
-template <typename D2=D2_Base&, typename D3=D3_Base&, typename DM=DM_Base&, typename C2=Cut_Base&, typename C3=Cut_Base&, typename CM=Cut_Base&>
+template <typename D2=tadah::models::D2_Base&, typename D3=tadah::models::D3_Base&, typename DM=tadah::models::DM_Base&, typename C2=tadah::models::Cut_Base&, typename C3=tadah::models::Cut_Base&, typename CM=tadah::models::Cut_Base&>
class DescriptorsCalc: public DC_Base {
public:
- Config &config;
+ tadah::core::Config &config;
/** Constructor to fully initialise this object.
*
* REQUIRED CONFIG KEYS:
@@ -61,7 +62,7 @@ class DescriptorsCalc: public DC_Base {
* at least one of: \ref RCUT2B, \ref RCUT3B,\ref RCUTMB
*
*/
- DescriptorsCalc(Config &c);
+ DescriptorsCalc(tadah::core::Config &c);
/** This constructor is equivalent to the main constructor above.
*
@@ -76,10 +77,10 @@ class DescriptorsCalc: public DC_Base {
*
* Usage example:
* \code{.cpp}
- * # Build calculator using Config object.
- * D2_Base *d2 = new D2_LJ(config);
- * D2_Base *d3 = new D3_Dummy(config);
- * D2_Base *dm = new DM_Dummy(config);
+ * # Build calculator using tadah::core::Config object.
+ * tadah::models::D2_Base *d2 = new D2_LJ(config);
+ * tadah::models::D2_Base *d3 = new D3_Dummy(config);
+ * tadah::models::D2_Base *dm = new DM_Dummy(config);
* DescriptorsCalc<> dc(config, d2, d3, dm);
*
* # Calculate all descriptors in a Structure st
@@ -101,16 +102,16 @@ class DescriptorsCalc: public DC_Base {
*
*/
template <typename T1, typename T2, typename T3,typename T4, typename T5, typename T6>
- DescriptorsCalc(Config &c, T1 &d2, T2 &d3, T3 &dm, T4 &c2, T5 &c3, T6 &cm);
+ DescriptorsCalc(tadah::core::Config &c, T1 &d2, T2 &d3, T3 &dm, T4 &c2, T5 &c3, T6 &cm);
/** Calculate all descriptors in a Structure st */
- StDescriptors calc(const Structure &st);
+ StDescriptors calc(const Structure &st, const StructureNeighborView &st_nb) override;
/** Calculate all descriptors in a StructureDB st_db */
- StDescriptorsDB calc(const StructureDB &st_db);
+ StDescriptorsDB calc(const StructureDB &st_db, const NeighborListDB &nldb) override;
/** Calculate density vector for the structure */
- void calc_rho(const Structure &st, StructureNeighborView &st_nb, StDescriptors &std);
+ void calc_rho(const Structure &st, const StructureNeighborView &st_nb, StDescriptors &std);
C2 c2;
C3 c3;
@@ -120,9 +121,9 @@ class DescriptorsCalc: public DC_Base {
DM dm;
private:
template <typename T1, typename T2, typename T3>
- DescriptorsCalc(Config &c,T1 &t1,T2 &t2,T3 &t3);
+ DescriptorsCalc(tadah::core::Config &c,T1 &t1,T2 &t2,T3 &t3);
- void calc(const Structure &st, StructureNeighborView &st_nb, StDescriptors &std);
+ void calc(const Structure &st, const StructureNeighborView &st_nb, StDescriptors &std);
void calc_dimer(const Structure &st, StDescriptors &std);
void common_constructor();
};
diff --git a/include/tadah/mlip/descriptors_calc.hpp b/include/tadah/mlip/descriptors_calc.hpp
index 1c47da2a3d599856a2563c627ce095506ad9a1ec..940ab70b060b66c5d53f8b6fe93b72e6ae6608b2 100644
--- a/include/tadah/mlip/descriptors_calc.hpp
+++ b/include/tadah/mlip/descriptors_calc.hpp
@@ -13,7 +13,7 @@ namespace tadah {
namespace mlip {
template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
template <typename T1, typename T2, typename T3>
-DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c, T1 &t1, T2 &t2, T3 &t3):
+DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(tadah::core::Config &c, T1 &t1, T2 &t2, T3 &t3):
config(c),
d2(t1),
d3(t2),
@@ -34,7 +34,7 @@ DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c, T1 &t1, T2 &t2, T
}
}
template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
-DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c):
+DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(tadah::core::Config &c):
DescriptorsCalc(c,c,c,c)
{
if (c.get<bool>("INIT2B")) {
@@ -80,7 +80,7 @@ DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c):
template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
template <typename T1, typename T2, typename T3,typename T4, typename T5, typename T6>
-DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(Config &c, T1 &d2, T2 &d3, T3 &dm, T4 &c2, T5 &c3, T6 &cm):
+DescriptorsCalc<D2,D3,DM,C2,C3,CM>::DescriptorsCalc(tadah::core::Config &c, T1 &d2, T2 &d3, T3 &dm, T4 &c2, T5 &c3, T6 &cm):
config(c),
c2(c2),
c3(c3),
@@ -164,13 +164,13 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::common_constructor() {
template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_rho(const Structure &st, const StructureNeighborView &st_nb, StDescriptors &st_d) {
double rcut_mb_sq = pow(config.get<double>("RCUTMBMAX"),2);
- rhos_type &rhos = st_d.rhos;
+ tadah::core::rhos_type &rhos = st_d.rhos;
size_t s = dm.rhoi_size()+dm.rhoip_size();
rhos.resize(s,st.natoms());
rhos.set_zero();
- Vec3d delij;
- Vec3d rj;
+ tadah::core::Vec3d delij;
+ tadah::core::Vec3d rj;
for (size_t i=0; i<st.natoms(); ++i) {
const Atom &a1 = st(i);
@@ -178,7 +178,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_rho(const Structure &st, const Str
const std::size_t *NbPtr = st_nb.getNeighborsPtr(i);
size_t NNbrs = st_nb.numNeighbors(i);
for (size_t jj=0; jj<NNbrs; ++jj) {
- // const Vec3d &a2pos = st.nn_pos(i,jj);
+ // const tadah::core::Vec3d &a2pos = st.nn_pos(i,jj);
size_t gj = NbPtr[jj];
size_t j = st.globalToLocal(gj); // j is local to st
@@ -228,7 +228,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const Structu
bias++;
if (initmb)
- calc_rho(st,st_d);
+ calc_rho(st,st_nb,st_d);
double rcut_max_sq = pow(config.get<double>("RCUTMAX"),2);
double rcut_2b_sq = 0.0;
@@ -239,7 +239,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const Structu
// zero all aeds and set bias
for (size_t i=0; i<st.natoms(); ++i) {
- aed_type &aed = st_d.get_aed(i);
+ tadah::core::aed_type &aed = st_d.get_aed(i);
aed.set_zero();
aed(0)=static_cast<double>(bias); // set bias
}
@@ -249,8 +249,8 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const Structu
// can be calculated
if (initmb) {
for (size_t i=0; i<st.natoms(); ++i) {
- aed_type &aed = st_d.get_aed(i);
- rho_type& rhoi = st_d.get_rho(i);
+ tadah::core::aed_type &aed = st_d.get_aed(i);
+ tadah::core::rho_type& rhoi = st_d.get_rho(i);
dm.calc_aed(rhoi,aed);
}
}
@@ -258,17 +258,17 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const Structu
bool use_force = config.get<bool>("FORCE");
bool use_stress = config.get<bool>("STRESS");
- Vec3d delij;
- Vec3d rj;
+ tadah::core::Vec3d delij;
+ tadah::core::Vec3d rj;
for (size_t i=0; i<st.natoms(); ++i) {
const Atom &a1 = st(i);
- aed_type &aed = st_d.get_aed(i);
+ tadah::core::aed_type &aed = st_d.get_aed(i);
int Zi = a1.atomicNumber();
const std::size_t *NbPtr = st_nb.getNeighborsPtr(i);
size_t NNbrs = st_nb.numNeighbors(i);
for (size_t jj=0; jj<NNbrs; ++jj) {
- // const Vec3d &a2pos = st.nn_pos(i,jj);
+ // const tadah::core::Vec3d &a2pos = st.nn_pos(i,jj);
size_t gj = NbPtr[jj];
size_t j = st.globalToLocal(gj); // j is local to st
@@ -296,7 +296,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const Structu
// CALCULATE TWO-BODY TERM
if (use_force || use_stress) {
- fd_type &fd_ij = st_d.fd[i][jj];
+ tadah::core::fd_type &fd_ij = st_d.fd[i][jj];
if (rij_sq <= rcut_2b_sq && init2b) {
d2.calc_all(Zi,Zj,rij,rij_sq,aed,fd_ij,0.5);
// Two-body descriptor calculates x-direction only - fd_ij(n,0)
@@ -311,7 +311,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const Structu
}
// CALCULATE MANY-BODY TERM
if (rij_sq <= rcut_mb_sq && initmb) {
- rho_type& rhoi = st_d.get_rho(i);
+ tadah::core::rho_type& rhoi = st_d.get_rho(i);
dm.calc_dXijdri(Zi,Zj,rij,rij_sq,delij,rhoi,fd_ij);
}
}
@@ -369,11 +369,11 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescr
// Not that this differ to how lammps implements this
// TODO make it consistent between those two
- Mat6R3C delM; //i1-j1,i2-j1,i1-j2,i2-j2,i1-i2,j1-j2
+ tadah::core::Mat6R3C delM; //i1-j1,i2-j1,i1-j2,i2-j2,i1-i2,j1-j2
double r_sq[6];
double r[6];
- const std::vector<Atom> &atoms = st.atoms;
+ /*const std::vector<Atom> &atoms = st.atoms;*/
size_t bias=0;
if (config.get<bool>("BIAS"))
@@ -396,7 +396,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescr
// zero all aeds+rho and set bias
for (size_t i=0; i<4; ++i) {
- aed_type &aed = st_d.get_aed(i);
+ tadah::core::aed_type &aed = st_d.get_aed(i);
aed.set_zero();
aed(0)=static_cast<double>(bias); // set bias
}
@@ -406,10 +406,12 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescr
st_d.rhos.set_zero();
}
- Vec3d xicom = 0.5*(atoms[0].position + atoms[1].position);
- Vec3d xjcom = 0.5*(atoms[2].position + atoms[3].position);
+ /*tadah::core::Vec3d xicom = 0.5*(atoms[0].position + atoms[1].position);*/
+ /*tadah::core::Vec3d xjcom = 0.5*(atoms[2].position + atoms[3].position);*/
+ tadah::core::Vec3d xicom = 0.5*(st.positionView(0) + st.positionView(1));
+ tadah::core::Vec3d xjcom = 0.5*(st.positionView(2) + st.positionView(3));
- Vec3d del_com = xicom-xjcom;;
+ tadah::core::Vec3d del_com = xicom-xjcom;;
double r_com_sq = del_com * del_com;
if (r_com_sq > rcut_com_sq) {
return;
@@ -417,7 +419,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescr
// compute all 6 distances
for (size_t n=0; n<6; ++n) {
- delM.row(n) = atoms[idx[n].first].position - atoms[idx[n].second].position;
+ delM.row(n) = st.positionView(idx[n].first) - st.positionView(idx[n].second);
r_sq[n] = delM.row(n) * delM.row(n);
r[n] = sqrt(r_sq[n]);
}
@@ -441,8 +443,8 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescr
// calculate many-body aed
if (initmb) {
for (size_t i=0; i<4; ++i) {
- aed_type &aed = st_d.get_aed(i);
- rho_type& rhoi = st_d.get_rho(i);
+ tadah::core::aed_type &aed = st_d.get_aed(i);
+ tadah::core::rho_type& rhoi = st_d.get_rho(i);
dm.calc_aed(rhoi,aed);
}
}
@@ -456,30 +458,39 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_dimer(const Structure &st, StDescr
//}
}
-template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
-StDescriptors DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st) {
+template <typename D2, typename D3, typename DM, typename C2, typename C3,
+ typename CM>
+StDescriptors DescriptorsCalc<D2, D3, DM, C2, C3, CM>::calc(
+ const Structure &st, const StructureNeighborView &st_nb) {
StDescriptors st_d(st, config);
- if (config.get<bool>("DIMER",0))
- calc_dimer(st,st_d);
- else
- calc(st,st_d);
+ if (config.get<bool>("DIMER", 0)) {
+ calc_dimer(st, st_d);
+ } else {
+ calc(st, st_nb, st_d);
+ }
return st_d;
}
-template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
-StDescriptorsDB DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const StructureDB &stdb) {
+template <typename D2, typename D3, typename DM, typename C2, typename C3,
+ typename CM>
+StDescriptorsDB
+DescriptorsCalc<D2, D3, DM, C2, C3, CM>::calc(const StructureDB &stdb, const NeighborListDB &nldb) {
StDescriptorsDB st_desc_db(stdb, config);
#ifdef _OPENMP
- #pragma omp parallel for
+#pragma omp parallel for
#endif
- for(size_t i=0; i<stdb.size(); ++i)
- if (config.get<bool>("DIMER",0))
- calc_dimer(stdb(i),st_desc_db(i));
- else
- calc(stdb(i),st_desc_db(i));
+ for (size_t i = 0; i < stdb.size(); ++i) {
+ StructureNeighborView st_nb(nldb.nlist(), stdb(i));
+ if (config.get<bool>("DIMER", 0)) {
+ calc_dimer(stdb(i), st_desc_db(i));
+ } else {
+ calc(stdb(i), st_nb, st_desc_db(i));
+ }
+ }
return st_desc_db;
}
+
}
}
#endif
diff --git a/include/tadah/mlip/descriptors_calc_base.h b/include/tadah/mlip/descriptors_calc_base.h
index d1b901a151f34a53273866c644bfcf23a281b9b5..e62bd76bfa3ac1abfc2132c8776f0e9d9b44de6e 100644
--- a/include/tadah/mlip/descriptors_calc_base.h
+++ b/include/tadah/mlip/descriptors_calc_base.h
@@ -5,14 +5,17 @@
#include <tadah/mlip/structure_db.h>
#include <tadah/mlip/st_descriptors.h>
#include <tadah/mlip/st_descriptors_db.h>
+#include <tadah/mlip/neighbor_list_db.h>
+#include <tadah/mlip/structure_neighbor_view.h>
+#include <tadah/mlip/neighbor_list_db.h>
namespace tadah {
namespace mlip {
class DC_Base {
public:
virtual ~DC_Base() {};
- virtual StDescriptors calc(const Structure &) { return StDescriptors(); };
- virtual StDescriptorsDB calc(const StructureDB &) {return StDescriptorsDB();};
+ virtual StDescriptors calc(const Structure &, const StructureNeighborView &) { return StDescriptors(); };
+ virtual StDescriptorsDB calc(const StructureDB &, const NeighborListDB &) {return StDescriptorsDB();};
};
}
}
diff --git a/include/tadah/mlip/design_matrix/design_matrix.h b/include/tadah/mlip/design_matrix/design_matrix.h
index c7dcaa79627a007d3ac197e0eadafc894c413f9d..4d5b7042120903582a63c8a30bc64d547c90f78a 100644
--- a/include/tadah/mlip/design_matrix/design_matrix.h
+++ b/include/tadah/mlip/design_matrix/design_matrix.h
@@ -1,7 +1,6 @@
#ifndef DESIGN_MATRIX_H
#define DESIGN_MATRIX_H
-#include "tadah/mlip/structure_neighbor_view.h"
#include <tadah/mlip/st_descriptors_db.h>
#include <tadah/mlip/structure_db.h>
#include <tadah/mlip/descriptors_calc.h>
@@ -21,7 +20,7 @@ namespace mlip {
class DesignMatrixBase {
public:
- static int phi_rows_num(Config &c, int nstruct_tot, int natoms_tot) {
+ static int phi_rows_num(tadah::core::Config &c, int nstruct_tot, int natoms_tot) {
bool force = c.get<bool>("FORCE");
bool stress = c.get<bool>("STRESS");
int rows = nstruct_tot;
@@ -32,7 +31,7 @@ public:
};
-/** Design Matrix \f$ \Phi \f$ and the corresponding target vector \f$ \mathrm{T} \f$
+/** Design tadah::core::Matrix \f$ \Phi \f$ and the corresponding target vector \f$ \mathrm{T} \f$
*
* An aggregation matrix composed of \phi_i matrices
*
@@ -68,7 +67,7 @@ public:
* - Stresses are scaled by 1/6.
* - Energies are scaled by 1/N_i
*
- * @tparam F Function_Base child -> BF_Base or Kern_Base
+ * @tparam F tadah::models::Function_Base child -> tadah::models::BF_Base or tadah::models::Kern_Base
*
*/
template <typename F>
@@ -86,17 +85,17 @@ public:
/** \brief This constructor fully initialise this object
*
- * This class is used to build a Design Matrix.
+ * This class is used to build a Design tadah::core::Matrix.
*
* Usage example:
* \code{.cpp}
- * Config config("Config");
- * BF_Linear bf(config);
+ * tadah::core::Config config("tadah::core::Config");
+ * tadah::models::BF_Linear bf(config);
* DesignMatrix<LinearKernel> desmat(bf, config);
*
* \endcode
*/
- DesignMatrix(F &f, Config &c, memory::IMLIPWorkspaceManager& workspaceManager)
+ DesignMatrix(F &f, tadah::core::Config &c, memory::IMLIPWorkspaceManager& workspaceManager)
: f(f),
workspaceManager_(&workspaceManager),
ownWorkspaceManager(false),
@@ -117,7 +116,7 @@ public:
}
// Constructor without workspaceManager_ parameter (delegating constructor)
- DesignMatrix(F &f, Config &c)
+ DesignMatrix(F &f, tadah::core::Config &c)
: DesignMatrix(f, c, *new memory::MLIPWorkspaceManager())
{
ownWorkspaceManager = true; // Set ownership flag
@@ -170,7 +169,7 @@ public:
/** \brief Calculate descriptors and build design matrix. */
template <typename DC>
- void build(const StructureDB &stdb, NeighbourListDB &nldb,
+ void build(const StructureDB &stdb, NeighborListDB &nldb,
Normaliser &norm,
DC &dc) {
@@ -244,7 +243,7 @@ public:
if (force) {
double fscale = 1;
if (scale) fscale = stdb(s).fweight*fweightglob/stdb(s).natoms()/3.0;///e_std_dev/f_std_dev;
- for (const Atom &a : stdb(s).atoms) {
+ for (const Atom &a : stdb(s)) {
ws->Tlabels(j) = 1;
ws->T(j++) = a.fx()*fscale;
ws->Tlabels(j) = 1;
@@ -269,14 +268,14 @@ public:
}
}
-phi_type &getPhi() { return ws->Phi; }
-t_type &getT() { return ws->T; }
-t_type &getTlabels() { return ws->Tlabels; }
+tadah::core::phi_type &getPhi() { return ws->Phi; }
+tadah::core::t_type &getT() { return ws->T; }
+tadah::core::t_type &getTlabels() { return ws->Tlabels; }
private:
memory::IMLIPWorkspaceManager* workspaceManager_ = nullptr;
bool ownWorkspaceManager = false;
- Config & config;
+ tadah::core::Config & config;
size_t rows = 0;
size_t cols = 0;
bool force;
@@ -307,8 +306,8 @@ private:
}
void compute_stdevs(const StructureDB &stdb) {
- t_type evec(stdb.size());
- Matrix svec(stdb.size(),6);
+ tadah::core::t_type evec(stdb.size());
+ tadah::core::Matrix svec(stdb.size(),6);
size_t num_forces=0;
for (size_t s=0; s<stdb.size(); ++s) {
evec(s) = stdb(s).energy/stdb(s).natoms();
@@ -337,12 +336,12 @@ private:
if (force) {
size_t j=0;
- t_type fvec(num_forces);
+ tadah::core::t_type fvec(num_forces);
for (size_t s=0; s<stdb.size(); ++s) {
- for (const Atom &a : stdb(s).atoms) {
- fvec(j++) = a.force(0);
- fvec(j++) = a.force(1);
- fvec(j++) = a.force(2);
+ for (const Atom &a : stdb(s)) {
+ fvec(j++) = a.fx();
+ fvec(j++) = a.fy();
+ fvec(j++) = a.fz();
}
}
//fvec /= e_std_dev;
diff --git a/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_base.h b/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_base.h
index 0cd2f6daccdf1c8ffb95dba6b240ae1647aa8e49..95b97cc7d36672029ca7099e985c55ee7830fca7 100644
--- a/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_base.h
+++ b/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_base.h
@@ -12,17 +12,17 @@
namespace tadah {
namespace mlip {
-struct DM_BF_Base: public DM_Function_Base, public virtual BF_Base {
+struct DM_BF_Base: public DM_Function_Base, public virtual tadah::models::BF_Base {
DM_BF_Base();
- DM_BF_Base(const Config &c);
+ DM_BF_Base(const tadah::core::Config &c);
virtual ~DM_BF_Base();
- // virtual size_t get_phi_cols(const Config &config)=0;
- // virtual void calc_phi_energy_row(phi_type &Phi, size_t &row,
+ // virtual size_t get_phi_cols(const tadah::core::Config &config)=0;
+ // virtual void calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row,
// const double fac, const Structure &st, const StDescriptors &st_d)=0;
- // virtual void calc_phi_force_rows(phi_type &Phi, size_t &row,
+ // virtual void calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row,
// const double fac, const Structure &st, const StDescriptors &st_d)=0;
- // virtual void calc_phi_stress_rows(phi_type &Phi, size_t &row,
+ // virtual void calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row,
// const double fac[6], const Structure &st, const StDescriptors &st_d)=0;
};
}
diff --git a/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_linear.h b/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_linear.h
index 96312bf6917cff054f7e1d64fb66b81c909d9101..a16c52ef2ae41eb2710da9bb717c46dc09b149c8 100644
--- a/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_linear.h
+++ b/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_linear.h
@@ -6,16 +6,16 @@
namespace tadah {
namespace mlip {
-struct DM_BF_Linear: public DM_BF_Base, public BF_Linear
+struct DM_BF_Linear: public DM_BF_Base, public tadah::models::BF_Linear
{
DM_BF_Linear();
- DM_BF_Linear(const Config &c);
- size_t get_phi_cols(const Config &config) override;
- void calc_phi_energy_row(phi_type &Phi, size_t &row,
+ DM_BF_Linear(const tadah::core::Config &c);
+ size_t get_phi_cols(const tadah::core::Config &config) override;
+ void calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- void calc_phi_force_rows(phi_type &Phi, size_t &row,
+ void calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- void calc_phi_stress_rows(phi_type &Phi, size_t &row,
+ void calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac[6], const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
};
}
diff --git a/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_polynomial2.h b/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_polynomial2.h
index 7bc796eed9f5ca4dc357e87de3127164e9aebeaf..3df7acd78583aeaa3bc780918eb75118407a67c3 100644
--- a/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_polynomial2.h
+++ b/include/tadah/mlip/design_matrix/functions/basis_functions/dm_bf_polynomial2.h
@@ -6,16 +6,16 @@
namespace tadah {
namespace mlip {
-struct DM_BF_Polynomial2: public DM_BF_Base, public BF_Polynomial2
+struct DM_BF_Polynomial2: public DM_BF_Base, public tadah::models::BF_Polynomial2
{
DM_BF_Polynomial2();
- DM_BF_Polynomial2(const Config &c);
- size_t get_phi_cols(const Config &config) override;
- void calc_phi_energy_row(phi_type &Phi, size_t &row,
+ DM_BF_Polynomial2(const tadah::core::Config &c);
+ size_t get_phi_cols(const tadah::core::Config &config) override;
+ void calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- void calc_phi_force_rows(phi_type &Phi, size_t &row,
+ void calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- void calc_phi_stress_rows(phi_type &Phi, size_t &row,
+ void calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac[6], const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
};
}
diff --git a/include/tadah/mlip/design_matrix/functions/dm_function_base.h b/include/tadah/mlip/design_matrix/functions/dm_function_base.h
index 8d40e9fab291e5421805ed39629663e684ebd78b..db4ecc2362cc2a0c754e1e7478cd2881181ff974 100644
--- a/include/tadah/mlip/design_matrix/functions/dm_function_base.h
+++ b/include/tadah/mlip/design_matrix/functions/dm_function_base.h
@@ -17,23 +17,23 @@
namespace tadah {
namespace mlip {
/** Base class for Kernels and Basis Functions */
-struct DM_Function_Base: public virtual Function_Base {
+struct DM_Function_Base: public virtual tadah::models::Function_Base {
- // Derived classes must implement Derived() and Derived(Config)
+ // Derived classes must implement Derived() and Derived(tadah::core::Config)
DM_Function_Base();
- DM_Function_Base(const Config &c);
+ DM_Function_Base(const tadah::core::Config &c);
virtual ~DM_Function_Base();
- virtual size_t get_phi_cols(const Config &)=0;
- virtual void calc_phi_energy_row(phi_type &, size_t &,
+ virtual size_t get_phi_cols(const tadah::core::Config &)=0;
+ virtual void calc_phi_energy_row(tadah::core::phi_type &, size_t &,
const double , const Structure &, const StructureNeighborView &st_nb, const StDescriptors &)=0;
- virtual void calc_phi_force_rows(phi_type &, size_t &,
+ virtual void calc_phi_force_rows(tadah::core::phi_type &, size_t &,
const double , const Structure &, const StructureNeighborView &st_nb, const StDescriptors &)=0;
- virtual void calc_phi_stress_rows(phi_type &, size_t &,
+ virtual void calc_phi_stress_rows(tadah::core::phi_type &, size_t &,
const double[6], const Structure &, const StructureNeighborView &st_nb, const StDescriptors &)=0;
};
}
}
-//template<> inline CONFIG::Registry<DM_Function_Base>::Map CONFIG::Registry<DM_Function_Base>::registry{};
-//template<> inline CONFIG::Registry<DM_Function_Base,Config&>::Map CONFIG::Registry<DM_Function_Base,Config&>::registry{};
+//template<> inline tadah::core::Registry<DM_Function_Base>::Map tadah::core::Registry<DM_Function_Base>::registry{};
+//template<> inline tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Map tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::registry{};
#endif
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_base.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_base.h
index 996dce6def190684c5559707f10e2ea602173313..a7ff86e519b326b1c1b98a3807d495de536cadab 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_base.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_base.h
@@ -20,16 +20,16 @@ namespace mlip {
* - ff = force descriptor
* - all derivatives are defined wrt to the second argument
*/
-class DM_Kern_Base: public DM_Function_Base, public virtual Kern_Base {
+class DM_Kern_Base: public DM_Function_Base, public virtual tadah::models::Kern_Base {
public:
DM_Kern_Base();
- DM_Kern_Base(const Config&c);
+ DM_Kern_Base(const tadah::core::Config&c);
virtual ~DM_Kern_Base();
- virtual size_t get_phi_cols(const Config &config) override;
- virtual void calc_phi_energy_row(phi_type &Phi, size_t &row, const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- virtual void calc_phi_force_rows(phi_type &Phi, size_t &row, const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- virtual void calc_phi_stress_rows(phi_type &Phi, size_t &row, const double fac[6], const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
+ virtual size_t get_phi_cols(const tadah::core::Config &config) override;
+ virtual void calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row, const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
+ virtual void calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row, const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
+ virtual void calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row, const double fac[6], const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
};
}
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_linear.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_linear.h
index 968478dae3e1b010c62c0cab87e63ef29f9501d4..f2bf10bf4ace6f7c2c198ea55c9e1dc9970305a7 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_linear.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_linear.h
@@ -14,19 +14,19 @@ namespace mlip {
* K(\mathbf{x}, \mathbf{y}) = \mathbf{x}^T \mathbf{y} = \sum_i x^{(i)} y^{(i)}
* \f]
*
- * @see Kern_Base BF_Linear
+ * @see tadah::models::Kern_Base tadah::models::BF_Linear
*/
-class DM_Kern_Linear : public DM_Kern_Base, public Kern_Linear {
+class DM_Kern_Linear : public DM_Kern_Base, public tadah::models::Kern_Linear {
public:
DM_Kern_Linear ();
- DM_Kern_Linear (const Config &c);
+ DM_Kern_Linear (const tadah::core::Config &c);
- size_t get_phi_cols(const Config &config) override;
- void calc_phi_energy_row(phi_type &Phi, size_t &row,
+ size_t get_phi_cols(const tadah::core::Config &config) override;
+ void calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- void calc_phi_force_rows(phi_type &Phi, size_t &row,
+ void calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
- void calc_phi_stress_rows(phi_type &Phi, size_t &row,
+ void calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac[6], const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) override;
};
}
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_lq.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_lq.h
index 704f25aac835da955af64453eae73b39f4378686..c6893bab383df1c727885a01b0f0510b8e4b708e 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_lq.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_lq.h
@@ -6,10 +6,10 @@
namespace tadah {
namespace mlip {
-class DM_Kern_LQ : public DM_Kern_Base, public Kern_LQ {
+class DM_Kern_LQ : public DM_Kern_Base, public tadah::models::Kern_LQ {
public:
DM_Kern_LQ();
- DM_Kern_LQ(const Config &c);
+ DM_Kern_LQ(const tadah::core::Config &c);
};
}
}
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_polynomial.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_polynomial.h
index ed8b4b47c7e078356356524f5f1c8527ddbc6386..bc3cf4f0fd664ec21bd25e0b2dfff7933c087755 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_polynomial.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_polynomial.h
@@ -6,10 +6,10 @@
namespace tadah {
namespace mlip {
-class DM_Kern_Polynomial : public DM_Kern_Base, public Kern_Polynomial {
+class DM_Kern_Polynomial : public DM_Kern_Base, public tadah::models::Kern_Polynomial {
public:
DM_Kern_Polynomial();
- DM_Kern_Polynomial(const Config &c);
+ DM_Kern_Polynomial(const tadah::core::Config &c);
};
}
}
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_quadratic.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_quadratic.h
index 5d977111c722cf30d29e2722fffb589a2c366d82..4b7de5df4558257140cc4e882787a9e8d0aecc1d 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_quadratic.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_quadratic.h
@@ -6,10 +6,10 @@
namespace tadah {
namespace mlip {
-class DM_Kern_Quadratic : public DM_Kern_Base, public Kern_Quadratic {
+class DM_Kern_Quadratic : public DM_Kern_Base, public tadah::models::Kern_Quadratic {
public:
DM_Kern_Quadratic();
- DM_Kern_Quadratic(const Config &c);
+ DM_Kern_Quadratic(const tadah::core::Config &c);
};
}
}
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_rbf.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_rbf.h
index d68c1141f652d865b2665c014890f99b55b540b2..a86b87e48ca155a820ae41da4886684d34fd64b7 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_rbf.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_rbf.h
@@ -6,10 +6,10 @@
namespace tadah {
namespace mlip {
-class DM_Kern_RBF : public DM_Kern_Base, public Kern_RBF {
+class DM_Kern_RBF : public DM_Kern_Base, public tadah::models::Kern_RBF {
public:
DM_Kern_RBF();
- DM_Kern_RBF(const Config &c);
+ DM_Kern_RBF(const tadah::core::Config &c);
};
}
}
diff --git a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_sigmoid.h b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_sigmoid.h
index 75c47a0b400bb85eb75cb581b513ce877db64278..ce0261ac8ae92c664d33fc7d5c82a9303f32293f 100644
--- a/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_sigmoid.h
+++ b/include/tadah/mlip/design_matrix/functions/kernels/dm_kern_sigmoid.h
@@ -6,10 +6,10 @@
namespace tadah {
namespace mlip {
-class DM_Kern_Sigmoid : public DM_Kern_Base, public Kern_Sigmoid {
+class DM_Kern_Sigmoid : public DM_Kern_Base, public tadah::models::Kern_Sigmoid {
public:
DM_Kern_Sigmoid();
- DM_Kern_Sigmoid(const Config &c);
+ DM_Kern_Sigmoid(const tadah::core::Config &c);
};
}
}
diff --git a/include/tadah/mlip/memory/DesignMatrixWorkspace.h b/include/tadah/mlip/memory/DesignMatrixWorkspace.h
index d9c9150e5ce86b8366b0d4458dc19a967b50e7f4..bf2d13ff51479e3d278553b70bd699c7dd3715c9 100644
--- a/include/tadah/mlip/memory/DesignMatrixWorkspace.h
+++ b/include/tadah/mlip/memory/DesignMatrixWorkspace.h
@@ -47,9 +47,9 @@ public:
*/
bool isSufficient(size_t m_, size_t n_) const;
- phi_type Phi;
- t_type T;
- t_type Tlabels; // 0-Energy, 1-Force, 2-Stress // TODO should be array of int not double
+ tadah::core::phi_type Phi;
+ tadah::core::t_type T;
+ tadah::core::t_type Tlabels; // 0-Energy, 1-Force, 2-Stress // TODO should be array of int not double
private:
diff --git a/include/tadah/mlip/models/basis.h b/include/tadah/mlip/models/basis.h
index d5b699b122276647e6a4fc5a686e158f9a72d1c8..73df373f37f33ae9a22503a69e8f76ad5aafbee8 100644
--- a/include/tadah/mlip/models/basis.h
+++ b/include/tadah/mlip/models/basis.h
@@ -17,18 +17,18 @@ namespace mlip {
template <typename K>
class Basis {
private:
- Config &config;
+ tadah::core::Config &config;
int verbose;
public:
- Matrix b;
- t_type T; // Vectors corresponding to basis vectors
- Basis(Config &c):
+ tadah::core::Matrix b;
+ tadah::core::t_type T; // Vectors corresponding to basis vectors
+ Basis(tadah::core::Config &c):
config(c),
verbose(c.get<int>("VERBOSE"))
{}
- void set_basis(Matrix &b_) {
+ void set_basis(tadah::core::Matrix &b_) {
b=b_;
}
void build_random_basis(size_t s, StDescriptorsDB &st_desc_db) {
@@ -59,7 +59,7 @@ larger than the amount of available AEDs\n");
for (size_t i=1; i<s; ++i) {
size_t st = std::get<0>(indices[i]);
size_t a = std::get<1>(indices[i]);
- const aed_type &aed = st_desc_db(st).get_aed(a);
+ const tadah::core::aed_type &aed = st_desc_db(st).get_aed(a);
for (size_t j=0; j<aed.size(); ++j) {
b(j,i)=aed[j];
}
@@ -100,7 +100,7 @@ larger than the amount of available AEDs\n");
const size_t st = indices[i];
T(i)=stdb(st).energy/st_desc_db(st).naed();
for( size_t a=0; a<st_desc_db(st).naed(); a++ ) {
- const aed_type &aed = st_desc_db(st).get_aed(a);
+ const tadah::core::aed_type &aed = st_desc_db(st).get_aed(a);
for (size_t j=0; j<aed.size(); ++j) {
b(j,i)+=aed[j]/st_desc_db(st).naed();
}
diff --git a/include/tadah/mlip/models/m_blr.h b/include/tadah/mlip/models/m_blr.h
index ceca9024b920f987d0e6c9211107e7576aa2df91..db4874b9a4562f00edf7549b1f656746a2b36247 100644
--- a/include/tadah/mlip/models/m_blr.h
+++ b/include/tadah/mlip/models/m_blr.h
@@ -41,27 +41,27 @@ namespace mlip {
*/
template
<class BF=DM_Function_Base&>
-class M_BLR: public M_Tadah_Base, public M_BLR_Train<BF> {
+class M_BLR: public M_Tadah_Base, public tadah::models::M_BLR_Train<BF> {
public:
- using M_BLR_Train<BF>::config;
- using M_BLR_Train<BF>::bf;
+ using tadah::models::M_BLR_Train<BF>::config;
+ using tadah::models::M_BLR_Train<BF>::bf;
/**
* @brief Initializes for training or prediction using a configuration.
*
* **Example**:
* \code{.cpp}
- * Config config("Config");
- * M_BLR<BF_Linear> blr(config);
+ * tadah::core::Config config("tadah::core::Config");
+ * M_BLR<tadah::models::BF_Linear> blr(config);
* \endcode
*
* @param c Configuration object.
*/
- M_BLR(Config &c):
- M_BLR_Train<BF>(c),
- desmat(M_BLR_Train<BF>::bf,c)
+ M_BLR(tadah::core::Config &c):
+ tadah::models::M_BLR_Train<BF>(c),
+ desmat(tadah::models::M_BLR_Train<BF>::bf,c)
{
norm = Normaliser(c);
}
@@ -72,29 +72,29 @@ public:
* @param bf Basis function.
* @param c Configuration object.
*/
- M_BLR(BF &bf, Config &c):
- M_BLR_Train<BF>(bf,c),
+ M_BLR(BF &bf, tadah::core::Config &c):
+ tadah::models::M_BLR_Train<BF>(bf,c),
desmat(bf,c)
{
norm = Normaliser(c);
}
- M_BLR(BF &bf, Config &c, tadah::mlip::memory::IMLIPWorkspaceManager& workspaceManager):
- M_BLR_Train<BF>(bf,c,workspaceManager),
+ M_BLR(BF &bf, tadah::core::Config &c, tadah::mlip::memory::IMLIPWorkspaceManager& workspaceManager):
+ tadah::models::M_BLR_Train<BF>(bf,c,workspaceManager),
desmat(bf,c,workspaceManager)
{
norm = Normaliser(c);
}
- double epredict(const aed_type &aed) const{
+ double epredict(const tadah::core::aed_type &aed) const{
return bf.epredict(weights,aed);
};
- double fpredict(const fd_type &fdij, const aed_type &aedi, const size_t k) const{
+ double fpredict(const tadah::core::fd_type &fdij, const tadah::core::aed_type &aedi, const size_t k) const{
return bf.fpredict(weights,fdij,aedi,k);
}
- force_type fpredict(const fd_type &fdij, const aed_type &aedi) const{
+ tadah::core::force_type fpredict(const tadah::core::fd_type &fdij, const tadah::core::aed_type &aedi) const{
return bf.fpredict(weights,fdij,aedi);
}
@@ -119,7 +119,7 @@ public:
config.add("FORCE", "false");
config.add("STRESS", "false");
- StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
+ StDescriptorsDB st_desc_db_temp = dc.calc(stdb, nldb);
if(config.template get<bool>("NORM")) {
norm = Normaliser(config);
@@ -137,18 +137,18 @@ public:
train(desmat);
}
- Structure predict(const Config &c, StDescriptors &std, const Structure &st) {
- if(config.template get<bool>("NORM") && !std.normalised && bf.get_label()!="BF_Linear")
+ Structure predict(const tadah::core::Config &c, StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb) {
+ if(config.template get<bool>("NORM") && !std.normalised && bf.get_label()!="tadah::models::BF_Linear")
norm.normalise(std);
- return M_Tadah_Base::predict(c,std,st);
+ return M_Tadah_Base::predict(c,std,st,st_nb);
}
- StructureDB predict(Config &c, const StructureDB &stdb, DC_Base &dc) {
- return M_Tadah_Base::predict(c,stdb,dc);
+ StructureDB predict(tadah::core::Config &c, const StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc) {
+ return M_Tadah_Base::predict(c,stdb,nldb,dc);
}
- Config get_param_file() {
- Config c = config;
+ tadah::core::Config get_param_file() {
+ tadah::core::Config c = config;
//c.remove("ALPHA");
//c.remove("BETA");
c.remove("DBFILE");
@@ -176,10 +176,10 @@ public:
c.clear_internal_keys();
return c;
}
- StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
- aed_type &predicted_error) {
+ StructureDB predict(tadah::core::Config config_pred, StructureDB &stdb, DC_Base &dc,
+ tadah::core::aed_type &predicted_error) {
- LinearRegressor::read_sigma(config_pred,Sigma);
+ tadah::models::LinearRegressor::read_sigma(config_pred,Sigma);
DesignMatrix<BF> dm(bf,config_pred);
dm.scale=false; // do not scale energy, forces and stresses
dm.build(stdb,norm,dc);
@@ -188,7 +188,7 @@ public:
double pmean = sqrt(predicted_error.mean());
// compute energy, forces and stresses
- aed_type Tpred = T_dgemv(dm.getPhi(), weights);
+ tadah::core::aed_type Tpred = T_dgemv(dm.getPhi(), weights);
// Construct StructureDB object with predicted values
StructureDB stdb_;
@@ -224,10 +224,10 @@ public:
if(!trained) throw std::runtime_error("This object is not trained!\n\
Hint: check different predict() methods.");
- phi_type &Phi = desmat.getPhi();
+ tadah::core::phi_type &Phi = desmat.getPhi();
// compute energy, forces and stresses
- aed_type Tpred = T_dgemv(Phi, weights);
+ tadah::core::aed_type Tpred = T_dgemv(Phi, weights);
double eweightglob=config.template get<double>("EWEIGHT");
double fweightglob=config.template get<double>("FWEIGHT");
@@ -269,12 +269,12 @@ private:
// normalise weights such that when predict is called
// we can supply it with a non-normalised descriptor
- t_type convert_to_nweights(const t_type &weights) const {
- if(bf.get_label()!="BF_Linear") {
+ tadah::core::t_type convert_to_nweights(const tadah::core::t_type &weights) const {
+ if(bf.get_label()!="tadah::models::BF_Linear") {
throw std::runtime_error("Cannot convert weights to nweights for\n\
non linear basis function\n");
}
- t_type nw(weights.rows());
+ tadah::core::t_type nw(weights.rows());
nw(0) = weights(0);
for (size_t i=1; i<weights.size(); ++i) {
@@ -289,13 +289,13 @@ non linear basis function\n");
return nw;
}
// The opposite of convert_to_nweights()
- t_type convert_to_weights(const t_type &nw) const {
- if(bf.get_label()!="BF_Linear") {
+ tadah::core::t_type convert_to_weights(const tadah::core::t_type &nw) const {
+ if(bf.get_label()!="tadah::models::BF_Linear") {
throw std::runtime_error("Cannot convert nweights to weights for\n\
non linear basis function\n");
}
// convert normalised weights back to "normal"
- t_type w(nw.rows());
+ tadah::core::t_type w(nw.rows());
w(0) = nw(0);
for (size_t i=1; i<nw.size(); ++i) {
if (norm.std_dev[i] > std::numeric_limits<double>::min())
@@ -314,25 +314,25 @@ non linear basis function\n");
// the train method destroys the Phi matrix
// In consequence, we cannot use it for quick prediction
// The simple solution, for now, is to make a copy of the Phi matrix
- //phi_type &Phi = desmat.Phi;
- phi_type Phi = desmat.getPhi();
- t_type T = desmat.getT();
- //t_type &T = desmat.T;
- M_BLR_Train<BF>::train(Phi,T);
+ //tadah::core::phi_type &Phi = desmat.Phi;
+ tadah::core::phi_type Phi = desmat.getPhi();
+ tadah::core::t_type T = desmat.getT();
+ //tadah::core::t_type &T = desmat.T;
+ tadah::models::M_BLR_Train<BF>::train(Phi,T);
if (config.template get<bool>("NORM") &&
- bf.get_label()=="BF_Linear") {
+ bf.get_label()=="tadah::models::BF_Linear") {
weights = convert_to_nweights(weights);
}
}
// Do we want to confuse user with those and make them public?
// Either way they must be stated as below to silence clang warning
- using M_BLR_Train<BF>::predict;
- using M_BLR_Train<BF>::train;
- using M_BLR_Train<BF>::trained;
- using M_BLR_Train<BF>::weights;
- using M_BLR_Train<BF>::Sigma;
+ using tadah::models::M_BLR_Train<BF>::predict;
+ using tadah::models::M_BLR_Train<BF>::train;
+ using tadah::models::M_BLR_Train<BF>::trained;
+ using tadah::models::M_BLR_Train<BF>::weights;
+ using tadah::models::M_BLR_Train<BF>::Sigma;
};
}
}
diff --git a/include/tadah/mlip/models/m_krr.h b/include/tadah/mlip/models/m_krr.h
index 609c3803a35613cb35c2f288ea6a30e6cb4c11ba..63f93b3792a2104d0f5741bc084227986f6265bf 100644
--- a/include/tadah/mlip/models/m_krr.h
+++ b/include/tadah/mlip/models/m_krr.h
@@ -44,7 +44,7 @@ namespace mlip {
template
<class K=DM_Function_Base&>
class M_KRR: public M_Tadah_Base,
- public M_KRR_Train<K>
+ public tadah::models::M_KRR_Train<K>
{
public:
@@ -54,14 +54,14 @@ public:
*
* **Example**:
* \code{.cpp}
- * Config config("Config");
- * M_KRR<Kern_Linear> krr(config);
+ * tadah::core::Config config("tadah::core::Config");
+ * M_KRR<tadah::models::Kern_Linear> krr(config);
* \endcode
*
* @param c Configuration object.
*/
- M_KRR(Config &c):
- M_KRR_Train<K>(c),
+ M_KRR(tadah::core::Config &c):
+ tadah::models::M_KRR_Train<K>(c),
basis(c),
desmat(kernel,c)
{
@@ -74,30 +74,30 @@ public:
* @param kernel Kernel function.
* @param c Configuration object.
*/
- M_KRR(K &kernel, Config &c):
- M_KRR_Train<K>(kernel,c),
+ M_KRR(K &kernel, tadah::core::Config &c):
+ tadah::models::M_KRR_Train<K>(kernel,c),
basis(c),
desmat(kernel,c)
{
norm = Normaliser(c);
}
- M_KRR(K &kernel, Config &c, tadah::mlip::memory::IMLIPWorkspaceManager& workspaceManager):
- M_KRR_Train<K>(kernel,c,workspaceManager),
+ M_KRR(K &kernel, tadah::core::Config &c, tadah::mlip::memory::IMLIPWorkspaceManager& workspaceManager):
+ tadah::models::M_KRR_Train<K>(kernel,c,workspaceManager),
basis(c),
desmat(kernel,c)
{
norm = Normaliser(c);
}
- double epredict(const aed_type &aed) const {
+ double epredict(const tadah::core::aed_type &aed) const {
return kernel.epredict(weights,aed);
};
- double fpredict(const fd_type &fdij, const aed_type &aedi, const size_t k) const {
+ double fpredict(const tadah::core::fd_type &fdij, const tadah::core::aed_type &aedi, const size_t k) const {
return kernel.fpredict(weights,fdij,aedi,k);
}
- force_type fpredict(const fd_type &fdij, const aed_type &aedi) const {
+ tadah::core::force_type fpredict(const tadah::core::fd_type &fdij, const tadah::core::aed_type &aedi) const {
return kernel.fpredict(weights,fdij,aedi);
}
@@ -130,9 +130,9 @@ public:
"This KRR implementation does exist: "\
+ std::to_string(modelN)+"\n");
}
- void train1(StructureDB &stdb, DC_Base &dc) {
+ void train1(StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc) {
// KRR implementation using EKM
- if(config.template get<bool>("NORM") || kernel.get_label()!="Kern_Linear") {
+ if(config.template get<bool>("NORM") || kernel.get_label()!="tadah::models::Kern_Linear") {
// either build basis or prep normaliser
std::string force=config.template get<std::string>("FORCE");
std::string stress=config.template get<std::string>("STRESS");
@@ -142,7 +142,7 @@ public:
config.add("FORCE", "false");
config.add("STRESS", "false");
- StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
+ StDescriptorsDB st_desc_db_temp = dc.calc(stdb,nldb);
if(config.template get<bool>("NORM")) {
norm = Normaliser(config);
@@ -155,7 +155,7 @@ public:
config.add("FORCE", force);
config.add("STRESS", stress);
- if (kernel.get_label()!="Kern_Linear") {
+ if (kernel.get_label()!="tadah::models::Kern_Linear") {
basis.build_random_basis(config.template get<size_t>("SBASIS"),st_desc_db_temp);
desmat.f.set_basis(basis.b);
kernel.set_basis(basis.b);
@@ -163,11 +163,11 @@ public:
ekm.configure(basis.b);
}
}
- desmat.build(stdb,norm,dc);
+ desmat.build(stdb,nldb,norm,dc);
train(desmat);
}
- void train2(StructureDB &stdb, DC_Base &dc) {
+ void train2(StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc) {
// NEW (BASIC) IMPLEMENTATION OF KRR //
std::string force=config.template get<std::string>("FORCE");
@@ -176,7 +176,7 @@ public:
config.remove("STRESS");
config.add("FORCE", "false");
config.add("STRESS", "false");
- StDescriptorsDB st_desc_db_temp = dc.calc(stdb);
+ StDescriptorsDB st_desc_db_temp = dc.calc(stdb, nldb);
if(config.template get<bool>("NORM")) {
norm = Normaliser(config);
norm.learn(st_desc_db_temp);
@@ -188,21 +188,21 @@ public:
config.add("STRESS", stress);
basis.prep_basis_for_krr(st_desc_db_temp,stdb);
kernel.set_basis(basis.b);
- M_KRR_Train<K>::train2(basis.b, basis.T);
+ tadah::models::M_KRR_Train<K>::train2(basis.b, basis.T);
}
- Structure predict(const Config &c, StDescriptors &std, const Structure &st) {
- if(config.template get<bool>("NORM") && !std.normalised && kernel.get_label()!="Kern_Linear")
+ Structure predict(const tadah::core::Config &c, StDescriptors &std, const Structure &st) {
+ if(config.template get<bool>("NORM") && !std.normalised && kernel.get_label()!="tadah::models::Kern_Linear")
norm.normalise(std);
return M_Tadah_Base::predict(c,std,st);
}
- StructureDB predict(Config &c, const StructureDB &stdb, DC_Base &dc) {
+ StructureDB predict(tadah::core::Config &c, const StructureDB &stdb, DC_Base &dc) {
return M_Tadah_Base::predict(c,stdb,dc);
}
- Config get_param_file() {
- Config c = config;
+ tadah::core::Config get_param_file() {
+ tadah::core::Config c = config;
c.remove("ALPHA");
c.remove("BETA");
c.remove("DBFILE");
@@ -229,7 +229,7 @@ public:
c.add("NSTDEV", norm.std_dev[i]);
}
}
- if (kernel.get_label()!="Kern_Linear") {
+ if (kernel.get_label()!="tadah::models::Kern_Linear") {
// dump basis to the config file file
// make sure keys are not accidently assigned
if (c.exist("SBASIS"))
@@ -246,10 +246,10 @@ public:
c.clear_internal_keys();
return c;
}
- StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
- aed_type &predicted_error) {
+ StructureDB predict(tadah::core::Config config_pred, StructureDB &stdb, DC_Base &dc,
+ tadah::core::aed_type &predicted_error) {
- LinearRegressor::read_sigma(config_pred,Sigma);
+ tadah::models::LinearRegressor::read_sigma(config_pred,Sigma);
DesignMatrix<K> dm(kernel,config_pred);
dm.scale=false; // do not scale energy, forces and stresses
dm.build(stdb,norm,dc);
@@ -259,7 +259,7 @@ public:
double pmean = sqrt(predicted_error.mean());
// compute energy, forces and stresses
- aed_type Tpred = T_dgemv(dm.getPhi(), weights);
+ tadah::core::aed_type Tpred = T_dgemv(dm.getPhi(), weights);
// Construct StructureDB object with predicted values
StructureDB stdb_;
@@ -295,10 +295,10 @@ public:
if(!trained) throw std::runtime_error("This object is not trained!\n\
Hint: check different predict() methods.");
- phi_type &Phi = desmat.getPhi();
+ tadah::core::phi_type &Phi = desmat.getPhi();
// compute energy, forces and stresses
- aed_type Tpred = T_dgemv(Phi, weights);
+ tadah::core::aed_type Tpred = T_dgemv(Phi, weights);
double eweightglob=config.template get<double>("EWEIGHT");
double fweightglob=config.template get<double>("FWEIGHT");
@@ -339,13 +339,13 @@ private:
Basis<K> basis;
DesignMatrix<K> desmat;
- t_type convert_to_nweights(const t_type &weights) const {
- if(kernel.get_label()!="Kern_Linear") {
+ tadah::core::t_type convert_to_nweights(const tadah::core::t_type &weights) const {
+ if(kernel.get_label()!="tadah::models::Kern_Linear") {
throw std::runtime_error("Cannot convert weights to nweights for\n\
non linear kernel\n");
}
- t_type kw(weights.rows());
- if(config.template get<bool>("NORM") && kernel.get_label()=="Kern_Linear") {
+ tadah::core::t_type kw(weights.rows());
+ if(config.template get<bool>("NORM") && kernel.get_label()=="tadah::models::Kern_Linear") {
// normalise weights such that when predict is called
// we can supply it with a non-normalised descriptor
kw.resize(weights.rows());
@@ -364,13 +364,13 @@ non linear kernel\n");
return kw;
}
// The opposite of convert_to_nweights()
- t_type convert_to_weights(const t_type &kw) const {
- if(kernel.get_label()!="Kern_Linear") {
+ tadah::core::t_type convert_to_weights(const tadah::core::t_type &kw) const {
+ if(kernel.get_label()!="tadah::models::Kern_Linear") {
throw std::runtime_error("Cannot convert nweights to weights for\n\
non linear kernel\n");
}
// convert normalised weights back to "normal"
- t_type w(kw.rows());
+ tadah::core::t_type w(kw.rows());
w(0) = kw(0);
for (size_t i=1; i<kw.size(); ++i) {
if (norm.std_dev[i] > std::numeric_limits<double>::min())
@@ -386,26 +386,26 @@ non linear kernel\n");
template <typename D>
void train(D &desmat) {
// TODO see comments in M_BLR
- phi_type Phi = desmat.getPhi();
- t_type T = desmat.getT();
- M_KRR_Train<K>::train(Phi,T);
+ tadah::core::phi_type Phi = desmat.getPhi();
+ tadah::core::t_type T = desmat.getT();
+ tadah::models::M_KRR_Train<K>::train(Phi,T);
if (config.template get<bool>("NORM") &&
- kernel.get_label()=="Kern_Linear") {
+ kernel.get_label()=="tadah::models::Kern_Linear") {
weights = convert_to_nweights(weights);
}
}
// Do we want to confuse user with those and make them public?
// Either way they must be stated as below to silence clang warning
- using M_KRR_Train<K>::predict;
- using M_KRR_Train<K>::train;
- using M_KRR_Train<K>::trained;
- using M_KRR_Train<K>::weights;
- using M_KRR_Train<K>::Sigma;
- using M_KRR_Train<K>::config;
- using M_KRR_Train<K>::kernel;
- using M_KRR_Train<K>::ekm;
+ using tadah::models::M_KRR_Train<K>::predict;
+ using tadah::models::M_KRR_Train<K>::train;
+ using tadah::models::M_KRR_Train<K>::trained;
+ using tadah::models::M_KRR_Train<K>::weights;
+ using tadah::models::M_KRR_Train<K>::Sigma;
+ using tadah::models::M_KRR_Train<K>::config;
+ using tadah::models::M_KRR_Train<K>::kernel;
+ using tadah::models::M_KRR_Train<K>::ekm;
};
}
}
diff --git a/include/tadah/mlip/models/m_tadah_base.h b/include/tadah/mlip/models/m_tadah_base.h
index 7ac57c8c5a46b0eed89810578baecd506a42b7e2..bad261c7231d3c33b9ca30a201b94bc630280127 100644
--- a/include/tadah/mlip/models/m_tadah_base.h
+++ b/include/tadah/mlip/models/m_tadah_base.h
@@ -17,8 +17,8 @@ namespace mlip {
/** This interface provides functionality required from all models.
*/
class M_Tadah_Base:
- public virtual M_Core,
- public virtual M_Predict
+ public virtual tadah::models::M_Core,
+ public virtual tadah::models::M_Predict
{
public:
@@ -27,58 +27,58 @@ public:
virtual ~M_Tadah_Base() {};
/** \brief Predict total energy for a set of atoms. */
- using M_Predict::epredict;
- using M_Predict::fpredict;
- using M_Core::predict;
+ using tadah::models::M_Predict::epredict;
+ using tadah::models::M_Predict::fpredict;
+ using tadah::models::M_Core::predict;
double epredict(const StDescriptors &std);
///** \brief Predict force between a pair of atoms in a k-direction. */
- //virtual double fpredict(const fd_type &fdij, const aed_type &aedi, size_t k)=0;
+ //virtual double fpredict(const tadah::core::fd_type &fdij, const tadah::core::aed_type &aedi, size_t k)=0;
///** \brief Predict force between a pair of atoms. */
- //virtual force_type fpredict(const fd_type &fdij,
- // const aed_type &aedi)=0;
+ //virtual tadah::core::force_type fpredict(const tadah::core::fd_type &fdij,
+ // const tadah::core::aed_type &aedi)=0;
/** \brief Predict total force on an atom a. */
- virtual void fpredict(const size_t a, force_type &v,
+ virtual void fpredict(const size_t a, tadah::core::force_type &v,
const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict total force on an atom a. */
- virtual force_type fpredict(const size_t a, const StDescriptors &std,
+ virtual tadah::core::force_type fpredict(const size_t a, const StDescriptors &std,
const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict energy, forces and stresses for the Structure st. */
- virtual Structure predict(const Config &c,
+ virtual Structure predict(const tadah::core::Config &c,
/*not const*/ StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict energy, forces and stresses for a set of Structures.
*
* Use precalculated descriptors.
*/
- virtual StructureDB predict(const Config &c,
+ virtual StructureDB predict(const tadah::core::Config &c,
/*not const*/ StDescriptorsDB &st_desc_db,
const StructureDB &stdb, const NeighborListDB &nldb);
/** \brief Predict energy, forces and stresses for a set of Structures. */
- virtual StructureDB predict(Config &c, const StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc);
+ virtual StructureDB predict(tadah::core::Config &c, const StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc);
/** \brief Predict virial stresses in a Structure st. */
- virtual stress_type spredict(const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
+ virtual tadah::core::stress_type spredict(const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict virial stresses in all Structures. */
- virtual stress_type spredict(const size_t a, const StDescriptors &std,
+ virtual tadah::core::stress_type spredict(const size_t a, const StDescriptors &std,
const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict both virial stresses and forces for a Structure */
virtual void stress_force_predict(const StDescriptors &std, Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict virial stress for atom a in a Structure st. */
- virtual void spredict(const size_t a, stress_type &s,
+ virtual void spredict(const size_t a, tadah::core::stress_type &s,
const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
/** Return potential file. */
- virtual Config get_param_file()=0;
+ virtual tadah::core::Config get_param_file()=0;
// TRAINING
@@ -99,13 +99,13 @@ public:
*/
virtual void train(StDescriptorsDB &, const StructureDB &) {};
- virtual StructureDB predict(Config config_pred, StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc,
- aed_type &predicted_error)=0;
+ virtual StructureDB predict(tadah::core::Config config_pred, StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc,
+ tadah::core::aed_type &predicted_error)=0;
virtual StructureDB predict(StructureDB &stdb, const NeighborListDB &nldb)=0;
};
}
}
-//template<> inline Registry<M_Tadah_Base,DM_Function_Base&,Config&>::Map Registry<M_Tadah_Base,DM_Function_Base&,Config&>::registry{};
+//template<> inline Registry<M_Tadah_Base,DM_Function_Base&,tadah::core::Config&>::Map Registry<M_Tadah_Base,DM_Function_Base&,tadah::core::Config&>::registry{};
#endif
diff --git a/include/tadah/mlip/neighbor_calc.h b/include/tadah/mlip/neighbor_calc.h
index a26b9a7d48df5b2c3213c5cf8293f77f39171e8f..a41d42d10e45b8e812c7f066d4e0609a32fe57dd 100644
--- a/include/tadah/mlip/neighbor_calc.h
+++ b/include/tadah/mlip/neighbor_calc.h
@@ -1,112 +1,644 @@
-#pragma once
-
-#include <tadah/mlip/neighbor_list.h>
-#include <vector>
-#include <cmath>
-
/**
- * @file neighbor_calc.h
- * @brief Demonstrates building a NeighborList using a naive O(N^2) approach.
+ * @file neighbor_calc.h
+ * @brief Implements a two-method neighbor calculation (naive and binned) analogous
+ * to the old NN_Finder, using the new StructureDB / NeighborList approach.
+ *
+ * The class provides:
+ * - A naive O(N^2) distance-check build for a single structure (coords-only).
+ * - A combined naive/binned build for an entire StructureDB, deciding per-structure
+ * whether to use bins (if cell dimensions ≥ cutoff) or the naive fallback otherwise.
+ *
+ * Binning uses a linked-cell approach, indexing each atom into one of the bins
+ * (nBinsA, nBinsB, nBinsC) derived from cell dimensions and the chosen cutoff.
+ * When any dimension is smaller than the cutoff, the method falls back to naive.
*
- * This example:
- * - Loops over all pairs (i,j).
- * - Checks if distance^2 < cutoff^2.
- * - If true, adds them as neighbors in both directions (i->j, j->i).
- * - Sets shift arrays to zero and uses mirrorIndex_ to store symmetrical offsets.
+ * This code relies on BLAS/LAPACK for matrix inversion when needed, and emphasizes
+ * contiguous data usage where possible for efficiency.
*/
+#pragma once
+#include <cmath>
+#include <algorithm>
+#include <stdexcept>
+#include <cstring>
+#include <chrono>
+#include <iostream>
+#include <vector>
+#include <numeric>
+#include <limits>
+
+#include <tadah/core/lapack.h> // If cblas_* or dgemm is needed
+#include <tadah/mlip/neighbor_list.h>
+#include <tadah/mlip/structure_db.h>
+
namespace tadah {
namespace mlip {
/**
* @class NeighborCalc
- * @brief Populates a NeighborList by scanning atom pairs for distances below cutoff.
+ * @brief Builds neighbor lists (per-structure or entire DB).
+ *
+ * This implementation supports:
+ * - A naive approach based on an O(N^2) distance check.
+ * - A binned approach using a linked-cell method for periodic boundaries.
+ * - A simple dimension check to switch between binned or naive for each structure.
+ *
+ * Example usage for a single structure (coords-only, no periodic shifts):
+ * @code
+ * NeighborList nlist;
+ * NeighborCalc calc;
+ * calc.build(nlist, coords, nAtoms, cutoff);
+ * @endcode
+ *
+ * Example usage for an entire structure DB:
+ * @code
+ * NeighborList nlist;
+ * StructureDB db(config);
+ * NeighborCalc calc;
+ * calc.build(db, nlist, cutoff); // sets up a combined HPC neighbor list for all structures
+ * @endcode
*/
class NeighborCalc
{
public:
/**
- * @brief Build the neighbor list with a naive O(N^2) distance check.
+ * @brief Naive O(N^2) build for a single, non-periodic set of coordinates.
+ *
+ * - Distance checked for all pairs (i < j).
+ * - If dist^2 < cutoff^2, the pair is added to both i's neighbor list and j's neighbor list.
+ * - Shifts are stored as zero because no cell data is provided.
+ * - mirrorIndex_ is populated to reference symmetrical pairs.
+ *
+ * Implementation details:
+ * 1) A dynamic adjacency buffer is used to collect neighbors for each atom.
+ * 2) This avoids reallocation or multiple passes to identify pair counts.
+ * 3) A second pass flattens adjacency into the SoA neighbor arrays (NeighborList).
*
* @param nlist The NeighborList to fill.
- * @param coords Pointer to x,y,z coords in [3*nAtoms].
+ * @param coords Pointer to x,y,z coords in [3 * nAtoms] layout.
* @param nAtoms Number of atoms in coords.
- * @param cutoff Distance cutoff for neighbor inclusion.
+ * @param cutoff Distance cutoff for inclusion in neighbor list.
*/
void build(NeighborList &nlist,
const double *coords,
std::size_t nAtoms,
double cutoff)
{
- std::vector<std::size_t> counts(nAtoms, 0);
-
- for (std::size_t i = 0; i < nAtoms; ++i) {
- for (std::size_t j = i + 1; j < nAtoms; ++j) {
- double dx = coords[3*j + 0] - coords[3*i + 0];
- double dy = coords[3*j + 1] - coords[3*i + 1];
- double dz = coords[3*j + 2] - coords[3*i + 2];
- double r2 = dx*dx + dy*dy + dz*dz;
- if (r2 < cutoff*cutoff) {
- counts[i]++;
- counts[j]++;
+ // Reason for adjacency: storing neighbors per atom for subsequent flattening
+ std::vector<std::vector<std::size_t>> adjacency(nAtoms);
+ double cutoffSq = cutoff * cutoff;
+
+ // Reserve adjacency to reduce reallocation
+ for (auto &vec : adjacency) {
+ vec.reserve(64);
+ }
+
+ // Pass1: collect neighbors in adjacency lists
+ // Only half-matrix is traversed (i < j) to avoid double counting
+ // Then symmetrical append is done immediately for j->i
+ for (std::size_t i = 0; i < nAtoms; i++) {
+ const double xi = coords[3*i + 0];
+ const double yi = coords[3*i + 1];
+ const double zi = coords[3*i + 2];
+ for (std::size_t j = i + 1; j < nAtoms; j++) {
+ double dx = xi - coords[3*j + 0];
+ double dy = yi - coords[3*j + 1];
+ double dz = zi - coords[3*j + 2];
+ double dist2 = dx*dx + dy*dy + dz*dz;
+ if(dist2 < cutoffSq) {
+ adjacency[i].push_back(j);
+ adjacency[j].push_back(i);
}
}
}
- std::size_t totalCount = 0;
- for (auto c : counts) {
- totalCount += c;
+ // Pass2: prefix sums to compute final neighbor-list memory usage
+ // prefixSum_[0] = 0, prefixSum_[i+1] = prefixSum_[i] + adjacency[i].size()
+ std::vector<std::size_t> prefixSum(nAtoms + 1, 0);
+ for (std::size_t i = 0; i < nAtoms; i++){
+ prefixSum[i+1] = prefixSum[i] + adjacency[i].size();
}
+ const std::size_t totalPairs = prefixSum[nAtoms];
- nlist.resize(nAtoms, totalCount);
+ // Resize HPC neighbor arrays
+ nlist.resize(nAtoms, totalPairs);
+ std::size_t *prefixPtr = nlist.prefixData();
+ std::size_t *neighPtr = nlist.neighborData();
+ int *sXPtr = nlist.shiftXData();
+ int *sYPtr = nlist.shiftYData();
+ int *sZPtr = nlist.shiftZData();
+ std::size_t *mirrorPtr = nlist.mirrorIndexData();
- std::size_t *pref = nlist.prefixData();
- pref[0] = 0;
- for (std::size_t i = 0; i < nAtoms; i++) {
- pref[i+1] = pref[i] + counts[i];
- }
-
- std::vector<std::size_t> writePos(nAtoms, 0);
-
- std::size_t *nbData = nlist.neighborData();
- int *sxData = nlist.shiftXData();
- int *syData = nlist.shiftYData();
- int *szData = nlist.shiftZData();
- std::size_t *mirrData = nlist.mirrorIndexData();
-
- for (std::size_t i = 0; i < nAtoms; ++i) {
- for (std::size_t j = i + 1; j < nAtoms; ++j) {
- double dx = coords[3*j + 0] - coords[3*i + 0];
- double dy = coords[3*j + 1] - coords[3*i + 1];
- double dz = coords[3*j + 2] - coords[3*i + 2];
- double r2 = dx*dx + dy*dy + dz*dz;
-
- if (r2 < cutoff*cutoff) {
- std::size_t offI = pref[i] + writePos[i];
- nbData[offI] = j;
- sxData[offI] = 0;
- syData[offI] = 0;
- szData[offI] = 0;
-
- std::size_t offJ = pref[j] + writePos[j];
- nbData[offJ] = i;
- sxData[offJ] = 0;
- syData[offJ] = 0;
- szData[offJ] = 0;
-
- mirrData[offI] = offJ;
- mirrData[offJ] = offI;
-
- writePos[i]++;
- writePos[j]++;
+ // Copy prefix sums
+ for (std::size_t i=0; i <= nAtoms; i++){
+ prefixPtr[i] = prefixSum[i];
+ }
+
+ // Pass3: fill SoA arrays (neighbors_, shift*, mirrorIndex_)
+ // The shifts are zero. mirrorIndex_ is determined by matching symmetrical pairs.
+ // This can be done in a single pass if a map is used, or done in two passes.
+ // A straightforward approach:
+ // - fill neighbor arrays
+ // - store partial mirror = (i) so that j->(i) can find the offset
+ // - after all are placed, do a second pass to fix mirrorIndex_ references
+ //
+ // Here, do: for each i: for each local neighbor "adj" => fill HPC
+ // Store a pointer offset for that neighbor. Then store i as a placeholder.
+
+ // Store partial data
+ for (std::size_t i = 0; i < nAtoms; i++){
+ const std::size_t base = prefixSum[i];
+ for (std::size_t k = 0; k < adjacency[i].size(); k++){
+ std::size_t j = adjacency[i][k];
+ neighPtr[ base + k ] = j;
+ sXPtr [ base + k ] = 0;
+ sYPtr [ base + k ] = 0;
+ sZPtr [ base + k ] = 0;
+ // Temporarily store "i" instead of the final mirror offset
+ mirrorPtr[ base + k ] = i;
+ }
+ }
+
+ // 2nd pass sets mirrorIndex_ properly.
+ // For the neighbor list of i, each neighbor is j.
+ // The mirror offset is the index in j's neighbor array that equals i.
+ // So find it by scanning adjacency[j] for i.
+ // Then the HPC offset = prefixSum[j] + that local index.
+ for (std::size_t i = 0; i < nAtoms; i++){
+ const std::size_t baseI = prefixSum[i];
+ for (std::size_t k = 0; k < adjacency[i].size(); k++){
+ std::size_t j = neighPtr[ baseI + k ];
+ std::size_t iInJ = 0;
+ // find i in adjacency[j]
+ // (Possible to store a map or do direct search for simplicity.)
+ const auto &adjJ = adjacency[j];
+ bool foundI = false;
+ for(std::size_t m = 0; m < adjJ.size(); m++){
+ if(adjJ[m] == i){
+ iInJ = m;
+ foundI = true;
+ break;
+ }
}
+ if(!foundI){
+ throw std::runtime_error("NeighborCalc::build - mirror index search error");
+ }
+ std::size_t baseJ = prefixSum[j];
+ // the local offset in j's neighbor array
+ std::size_t mirrorOffset = baseJ + iInJ;
+ mirrorPtr[ baseI + k ] = mirrorOffset;
}
}
}
- void build (const StructureDB &stdb, NeighborList &nList, double cutoff) {
- }
+
+ /**
+ * @brief Builds a neighbor list for an entire StructureDB, deciding for each structure
+ * whether to do a binned or naive approach, then concatenating them in HPC form.
+ *
+ * - Each structure has an SoA offset sOff = sRef.offset().
+ * - The total HPC arrays span [0..(sum of all atoms in the DB)-1].
+ * - Bins are built individually within each structure. Inter-structure neighbors are not added.
+ * - If a structure’s cell dimension < cutoff, a naive fallback is used.
+ * - The final neighbor list captures all structures’ adjacency in a single SoA.
+ *
+ * Example:
+ * @code
+ * StructureDB db(config);
+ * NeighborList nlist;
+ * NeighborCalc calc;
+ * calc.build(db, nlist, 5.0);
+ * @endcode
+ */
+ void build(const StructureDB &stdb, NeighborList &nList, double cutoff)
+ {
+ // Identify total number of atoms
+ const std::size_t Nall = stdb.X_.size();
+ // Prepare an adjacency buffer for all atoms
+ std::vector<std::vector<std::size_t>> adjacency(Nall);
+ std::vector<std::vector<int>> shiftX(Nall), shiftY(Nall), shiftZ(Nall);
+
+ // Reserve typical neighbor capacity
+#pragma omp parallel for
+ for (std::size_t i = 0; i < Nall; i++){
+ adjacency[i].reserve(64);
+ shiftX[i].reserve(64);
+ shiftY[i].reserve(64);
+ shiftZ[i].reserve(64);
+ }
+
+ double cutoffSq = cutoff * cutoff;
+
+ // Loop over structures, each structure is processed individually,
+ // The HPC adjacency is appended at the global offsets specific to that structure.
+ // This forms a blockwise approach: no cross-structure neighbors.
+ std::size_t nStructs = stdb.size();
+
+#ifdef _OPENMP
+ // Potential to parallelize over structures if desired
+#pragma omp parallel for
+#endif
+ for (std::size_t sIndex = 0; sIndex < nStructs; sIndex++){
+ const Structure &st = stdb(sIndex);
+ // st.offset() is the global start index for this structure
+ const std::size_t sOff = st.offset();
+ const std::size_t sNatom= st.natoms();
+
+ // Decide binned or naive
+ if(!structureBoxIsLargeEnough(st, cutoffSq)) {
+ buildNaiveForStructure(st, adjacency, shiftX, shiftY, shiftZ, cutoffSq);
+ } else {
+ buildBinnedForStructure(st, adjacency, shiftX, shiftY, shiftZ, cutoffSq);
+ }
+ }
+
+ // Now adjacency is filled for the entire DB, structure by structure.
+ // Next step: prefix sums
+ std::vector<std::size_t> prefixSum(Nall + 1, 0);
+ for (std::size_t i = 0; i < Nall; i++){
+ prefixSum[i+1] = prefixSum[i] + adjacency[i].size();
+ }
+ const std::size_t totalPairs = prefixSum[Nall];
+
+ // Resize HPC arrays
+ nList.resize(Nall, totalPairs);
+ std::size_t *prefixPtr = nList.prefixData();
+ std::size_t *neighPtr = nList.neighborData();
+ int *sXPtr = nList.shiftXData();
+ int *sYPtr = nList.shiftYData();
+ int *sZPtr = nList.shiftZData();
+ std::size_t *mirrorPtr = nList.mirrorIndexData();
+
+ // Copy prefix sums
+ for (std::size_t i=0; i <= Nall; i++){
+ prefixPtr[i] = prefixSum[i];
+ }
+
+ // Fill HPC arrays (neighbors, shifts, mirror placeholders)
+#pragma omp parallel for
+ for (std::size_t i = 0; i < Nall; i++){
+ const std::size_t base = prefixSum[i];
+ const std::size_t deg = adjacency[i].size();
+ for(std::size_t k = 0; k < deg; k++){
+ neighPtr [base + k] = adjacency[i][k];
+ sXPtr [base + k] = shiftX[i][k];
+ sYPtr [base + k] = shiftY[i][k];
+ sZPtr [base + k] = shiftZ[i][k];
+ // Temporarily store "i" as placeholder
+ mirrorPtr[base + k] = i;
+ }
+ }
+
+ // Final pass to compute mirror offsets
+ // For neighbor i->j, the mirror offset is the index of i in adjacency[j].
+#pragma omp parallel for
+ for (std::size_t i = 0; i < Nall; i++){
+ const std::size_t baseI = prefixSum[i];
+ const std::size_t degI = adjacency[i].size();
+ for (std::size_t k = 0; k < degI; k++){
+ std::size_t j = neighPtr[baseI + k];
+ // Search adjacency[j] for i
+ const auto &adjJ = adjacency[j];
+ std::size_t foundIndex = 0;
+ bool found = false;
+ for (std::size_t m = 0; m < adjJ.size(); m++){
+ if (adjJ[m] == i){
+ foundIndex = m;
+ found = true;
+ break;
+ }
+ }
+ if(!found){
+ throw std::runtime_error("NeighborCalc::build - mismatch in adjacency mirror pass");
+ }
+ std::size_t jBase = prefixSum[j];
+ std::size_t mirrorOffs = jBase + foundIndex;
+ mirrorPtr[baseI + k] = mirrorOffs;
+ }
+ }
+ }
+
+private:
+ /**
+ * @brief Checks if all cell basis vectors are of length^2 >= cutoff^2.
+ * If any dimension < cutoff, returns false => naive approach recommended.
+ *
+ * @param st The structure in question.
+ * @param cutoffSq The squared cutoff for comparison.
+ */
+ bool structureBoxIsLargeEnough(const Structure &st, double cutoffSq) const
+ {
+ // Access each row of st.cell
+ // Reason: if row's length^2 < cutoff^2 => fallback to naive
+ for (int i = 0; i < 3; i++) {
+ double vx = st.cell(i, 0);
+ double vy = st.cell(i, 1);
+ double vz = st.cell(i, 2);
+ double len2 = vx*vx + vy*vy + vz*vz;
+ if (len2 < cutoffSq) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ /**
+ * @brief Naive approach for a single structure within the DB:
+ * - loops over (a1 < a2) in [sOff.. sOff+sNatom),
+ * - accounts for periodic images by enumerating possible (n1,n2,n3) if needed,
+ * - but here, it is coherent with the old NNFinder "calc_naive".
+ *
+ * Implementation detail:
+ * - For each pair, check distance including shift if (n1,n2,n3) != (0,0,0).
+ * - If distance^2 < cutoff^2, add adjacency in both directions.
+ * - shift arrays store integer (n1,n2,n3).
+ */
+ void buildNaiveForStructure(const Structure &st,
+ std::vector<std::vector<std::size_t>> &adjacency,
+ std::vector<std::vector<int>> &shiftX,
+ std::vector<std::vector<int>> &shiftY,
+ std::vector<std::vector<int>> &shiftZ,
+ double cutoffSq) const
+ {
+ const std::size_t sOff = st.offset();
+ const std::size_t sNatom = st.natoms();
+
+ // Determine how many images in ± direction to consider
+ // This replicates the logic from num_shifts in the old code.
+ int N[3];
+ computeNumShifts(st, N, std::sqrt(cutoffSq));
+
+ // Precompute real-space displacement of each shift
+ // plus the integer (n1,n2,n3) itself
+ std::vector<tadah::core::Vec3d> shifts;
+ std::vector<tadah::core::Vec3i> shiftIdx;
+ {
+ const int Nx = 2*N[0] + 1;
+ const int Ny = 2*N[1] + 1;
+ const int Nz = 2*N[2] + 1;
+ shifts.reserve(Nx * Ny * Nz);
+ shiftIdx.reserve(Nx * Ny * Nz);
+ for (int n1 = -N[0]; n1 <= N[0]; n1++){
+ for (int n2 = -N[1]; n2 <= N[1]; n2++){
+ for (int n3 = -N[2]; n3 <= N[2]; n3++){
+ tadah::core::Vec3d disp;
+ disp[0] = st.cell(0,0)*n1 + st.cell(1,0)*n2 + st.cell(2,0)*n3;
+ disp[1] = st.cell(0,1)*n1 + st.cell(1,1)*n2 + st.cell(2,1)*n3;
+ disp[2] = st.cell(0,2)*n1 + st.cell(1,2)*n2 + st.cell(2,2)*n3;
+ shifts.push_back(disp);
+ shiftIdx.push_back(tadah::core::Vec3i(n1, n2, n3));
+ }
+ }
+ }
+ }
+
+ // Shortcut references to SoA
+ const double *X = st.db_->X_.data();
+ const double *Y = st.db_->Y_.data();
+ const double *Z = st.db_->Z_.data();
+
+ // For each shift (except the self-shift, where a2 starts at a1+1),
+ // check distance^2 < cutoff^2 => store adjacency
+ for (std::size_t s = 0; s < shifts.size(); s++){
+ const tadah::core::Vec3d &sh = shifts[s];
+ const tadah::core::Vec3i &shI= shiftIdx[s];
+
+ bool isSelfShift = (shI[0]==0 && shI[1]==0 && shI[2]==0);
+ std::size_t startA2 = isSelfShift ? 1ul : 0ul;
+
+ for (std::size_t locA1 = 0; locA1 < sNatom; locA1++){
+ std::size_t gA1 = sOff + locA1; // global index
+ double x1 = X[gA1];
+ double y1 = Y[gA1];
+ double z1 = Z[gA1];
+ for (std::size_t locA2 = locA1 + startA2; locA2 < sNatom; locA2++){
+ std::size_t gA2 = sOff + locA2;
+ double dx = x1 - (X[gA2] + sh[0]);
+ double dy = y1 - (Y[gA2] + sh[1]);
+ double dz = z1 - (Z[gA2] + sh[2]);
+ double dist2 = dx*dx + dy*dy + dz*dz;
+ if(dist2 < cutoffSq){
+ adjacency[gA1].push_back(gA2);
+ shiftX[gA1].push_back(shI[0]);
+ shiftY[gA1].push_back(shI[1]);
+ shiftZ[gA1].push_back(shI[2]);
+
+ adjacency[gA2].push_back(gA1);
+ shiftX[gA2].push_back(-shI[0]);
+ shiftY[gA2].push_back(-shI[1]);
+ shiftZ[gA2].push_back(-shI[2]);
+ }
+ }
+ }
+ }
+ }
+
+ /**
+ * @brief Binning-based approach for a single structure.
+ * - Divides the cell into (nBinsA,nBinsB,nBinsC).
+ * - Each atom is placed in a bin, using fractional coordinates mod 1.
+ * - Then searches 3x3x3 adjacent bin cells for neighbors.
+ * - If any dimension < cutoff, naive fallback is indicated from structureBoxIsLargeEnough.
+ */
+ void buildBinnedForStructure(const Structure &st,
+ std::vector<std::vector<std::size_t>> &adjacency,
+ std::vector<std::vector<int>> &shiftX,
+ std::vector<std::vector<int>> &shiftY,
+ std::vector<std::vector<int>> &shiftZ,
+ double cutoffSq) const
+ {
+ const std::size_t sOff = st.offset();
+ const std::size_t sNatom = st.natoms();
+
+ // Invert cell
+ double invCell[9];
+ invertCell3x3(st.cell.data(), invCell);
+
+ // Compute number of bins per dimension
+ // nBinsX = floor( length(cellRow0)/cutoff ), etc. with at least 1 bin
+ auto rowLength = [&](int row){
+ double vx = st.cell(row, 0);
+ double vy = st.cell(row, 1);
+ double vz = st.cell(row, 2);
+ return std::sqrt(vx*vx + vy*vy + vz*vz);
+ };
+ double cellLenA = rowLength(0);
+ double cellLenB = rowLength(1);
+ double cellLenC = rowLength(2);
+
+ int nBinsA = std::max(1, int(std::floor(cellLenA / std::sqrt(cutoffSq))));
+ int nBinsB = std::max(1, int(std::floor(cellLenB / std::sqrt(cutoffSq))));
+ int nBinsC = std::max(1, int(std::floor(cellLenC / std::sqrt(cutoffSq))));
+
+ // Prepare bins
+ struct BinCell {
+ std::vector<std::size_t> atomIndices;
+ };
+ std::vector<BinCell> bins(nBinsA * nBinsB * nBinsC);
+
+ auto wrapBinIndex = [&](int x, int nbins) {
+ // Reason: ensures periodic wrap for negative or large indices
+ return ( (x % nbins) + nbins ) % nbins;
+ };
+ auto binIndex = [&](int ia, int ib, int ic){
+ ia = wrapBinIndex(ia, nBinsA);
+ ib = wrapBinIndex(ib, nBinsB);
+ ic = wrapBinIndex(ic, nBinsC);
+ return std::size_t(ia + nBinsA*(ib + nBinsB*ic));
+ };
+
+ // Fill bins with atoms
+ const double *X = st.db_->X_.data();
+ const double *Y = st.db_->Y_.data();
+ const double *Z = st.db_->Z_.data();
+
+ for (std::size_t locA = 0; locA < sNatom; locA++){
+ std::size_t gA = sOff + locA;
+ double px = X[gA];
+ double py = Y[gA];
+ double pz = Z[gA];
+ // fractional coords
+ double fx = invCell[0]*px + invCell[3]*py + invCell[6]*pz;
+ double fy = invCell[1]*px + invCell[4]*py + invCell[7]*pz;
+ double fz = invCell[2]*px + invCell[5]*py + invCell[8]*pz;
+ fx -= std::floor(fx);
+ fy -= std::floor(fy);
+ fz -= std::floor(fz);
+
+ int ia = int(std::floor(fx * nBinsA));
+ int ib = int(std::floor(fy * nBinsB));
+ int ic = int(std::floor(fz * nBinsC));
+
+ bins[ binIndex(ia, ib, ic) ].atomIndices.push_back(gA);
+ }
+
+ // Search neighbors among bin pairs in a 3x3x3 region
+ for (int ia = 0; ia < nBinsA; ia++){
+ for (int ib = 0; ib < nBinsB; ib++){
+ for (int ic = 0; ic < nBinsC; ic++){
+ std::size_t b0 = binIndex(ia, ib, ic);
+ auto &vec0 = bins[b0].atomIndices;
+ for (int ja = ia-1; ja <= ia+1; ja++){
+ for (int jb = ib-1; jb <= ib+1; jb++){
+ for (int jc = ic-1; jc <= ic+1; jc++){
+ std::size_t b1 = binIndex(ja, jb, jc);
+ auto &vec1 = bins[b1].atomIndices;
+ int dA = ja - ia;
+ int dB = jb - ib;
+ int dC = jc - ic;
+
+ // For each pair in these bins
+ for (auto gA : vec0){
+ for (auto gB : vec1){
+ // Avoid double counting
+ // If b0 == b1, only handle gB > gA
+ if(b0 == b1 && gB <= gA) {
+ continue;
+ }
+ double dx = X[gA] - (X[gB] + dA*(st.cell(0,0)/double(nBinsA))
+ + dB*(st.cell(0,1)/double(nBinsB))
+ + dC*(st.cell(0,2)/double(nBinsC)));
+ double dy = Y[gA] - (Y[gB] + dA*(st.cell(1,0)/double(nBinsA))
+ + dB*(st.cell(1,1)/double(nBinsB))
+ + dC*(st.cell(1,2)/double(nBinsC)));
+ double dz = Z[gA] - (Z[gB] + dA*(st.cell(2,0)/double(nBinsA))
+ + dB*(st.cell(2,1)/double(nBinsB))
+ + dC*(st.cell(2,2)/double(nBinsC)));
+ double r2 = dx*dx + dy*dy + dz*dz;
+ if(r2 < cutoffSq) {
+ adjacency[gA].push_back(gB);
+ shiftX[gA].push_back(dA);
+ shiftY[gA].push_back(dB);
+ shiftZ[gA].push_back(dC);
+
+ adjacency[gB].push_back(gA);
+ shiftX[gB].push_back(-dA);
+ shiftY[gB].push_back(-dB);
+ shiftZ[gB].push_back(-dC);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /**
+ * @brief Inverts a 3×3 cell matrix using LAPACK or direct cofactor approach (for speed).
+ * The input and output arrays are 3×3 in column-major order.
+ */
+ void invertCell3x3(const double *inMat, double *outMat) const
+ {
+ // If using an external LAPACK routine, e.g. dgetrf + dgetri, place code here.
+ // Alternatively, replicate the direct cofactor approach for speed:
+ // This internal method is adapted from NN_Finder's inverse_3x3_direct.
+ // Passive voice and personal pronouns avoided, as requested.
+ double tmp[9];
+ // Unrolled, but for clarity, gather into a local buffer
+ for(int i=0; i<9; i++){
+ tmp[i] = inMat[i];
+ }
+
+ // Cofactor-based inversion (same as old code):
+ // The numeric threshold can be tuned.
+ double m00 = tmp[0]; double m10 = tmp[3]; double m20 = tmp[6];
+ double m01 = tmp[1]; double m11 = tmp[4]; double m21 = tmp[7];
+ double m02 = tmp[2]; double m12 = tmp[5]; double m22 = tmp[8];
+ double det = m00*(m11*m22 - m12*m21)
+ - m01*(m10*m22 - m12*m20)
+ + m02*(m10*m21 - m11*m20);
+ if(std::fabs(det) < 1e-14) {
+ throw std::runtime_error("NeighborCalc::invertCell3x3: near-singular matrix");
+ }
+ double invDet = 1.0/det;
+ double i00 = (m11*m22 - m12*m21) * invDet;
+ double i01 = -(m01*m22 - m02*m21) * invDet;
+ double i02 = (m01*m12 - m02*m11) * invDet;
+ double i10 = -(m10*m22 - m12*m20) * invDet;
+ double i11 = (m00*m22 - m02*m20) * invDet;
+ double i12 = -(m00*m12 - m02*m10) * invDet;
+ double i20 = (m10*m21 - m11*m20) * invDet;
+ double i21 = -(m00*m21 - m01*m20) * invDet;
+ double i22 = (m00*m11 - m01*m10) * invDet;
+
+ outMat[0]=i00; outMat[1]=i01; outMat[2]=i02;
+ outMat[3]=i10; outMat[4]=i11; outMat[5]=i12;
+ outMat[6]=i20; outMat[7]=i21; outMat[8]=i22;
+ }
+
+ /**
+ * @brief Determines the integer image shifts needed for naive approach,
+ * as in old code's num_shifts().
+ *
+ * - Inverts cell to get reciprocal lengths l1, l2, l3
+ * - The shift bound N[i] is (cutoff * b_i / f_i),
+ * where b_i = max(1, floor( 1 / l_i )), f_i ~ 1/l_i
+ */
+ void computeNumShifts(const Structure &st, int N[3], double cutoff) const
+ {
+ // Invert st.cell using built-in function or direct approach
+ tadah::core::Matrix3d invCell = st.cell.inverse();
+ double l1 = invCell.col(0).norm();
+ double l2 = invCell.col(1).norm();
+ double l3 = invCell.col(2).norm();
+
+ double f1 = (l1 > 0.0) ? 1.0/l1 : 1.0;
+ double f2 = (l2 > 0.0) ? 1.0/l2 : 1.0;
+ double f3 = (l3 > 0.0) ? 1.0/l3 : 1.0;
+
+ // b_i = max(1, int(f_i / cutoff)) in old code
+ // Then N[i] ~ round( cutoff * b_i / f_i )
+ int b1 = std::max(int(f1/cutoff), 1);
+ int b2 = std::max(int(f2/cutoff), 1);
+ int b3 = std::max(int(f3/cutoff), 1);
+
+ N[0] = (int)(std::round(0.5 + cutoff*b1/f1));
+ N[1] = (int)(std::round(0.5 + cutoff*b2/f2));
+ N[2] = (int)(std::round(0.5 + cutoff*b3/f3));
+ }
};
} // end namespace mlip
} // end namespace tadah
-
diff --git a/include/tadah/mlip/nn_finder.h b/include/tadah/mlip/nn_finder.h
deleted file mode 100644
index 561d267fb6fd564615b901181386390d74f5878f..0000000000000000000000000000000000000000
--- a/include/tadah/mlip/nn_finder.h
+++ /dev/null
@@ -1,166 +0,0 @@
-#ifndef NN_FINDER_H
-#define NN_FINDER_H
-
-#include <tadah/core/config.h>
-#include <tadah/core/lapack.h>
-#include <tadah/mlip/structure.h>
-#include <tadah/mlip/structure_db.h>
-
-namespace tadah {
-namespace mlip {
-/**
- * @class NNFinder
- *
- * Nearest Neighbor Finder that constructs a full nearest neighbor list
- * for every atom in a structure, using:
- * - a **binned** (linked-cell) approach if every cell dimension ≥ cutoff,
- * - a **naive** approach otherwise (fallback).
- *
- * The cutoff is from config("RCUTMAX").
- */
-class NNFinder {
-private:
- double cutoff; ///< The cutoff distance
- double cutoff_sq; ///< cutoff^2 for distance checks
-
- /**
- * Return false if any cell dimension < cutoff => fallback to naive.
- */
- bool check_box(Structure &st);
-
- /**
- * Invert a 3×3 matrix (column-major)
- */
-/**
- * Invert a 3x3 matrix using the cofactor (adjugate) method,
- * matching the layout:
- * M(r,c) stored at index [r*3 + c]
- *
- * So, for example:
- * M(0,0) -> inM[0],
- * M(1,0) -> inM[3],
- * M(2,0) -> inM[6],
- * M(0,1) -> inM[1],
- * ...
- *
- * The output is also stored in the same pattern.
- *
- * Throws std::runtime_error if the matrix is near-singular.
- */
-inline void inverse_3x3_direct(const double* inM, double* outM)
-{
- // Map input array -> matrix elements:
- // M(r,c) => inM[r*3 + c]
- const double m00 = inM[0]; // M(0,0)
- const double m10 = inM[3]; // M(1,0)
- const double m20 = inM[6]; // M(2,0)
-
- const double m01 = inM[1]; // M(0,1)
- const double m11 = inM[4]; // M(1,1)
- const double m21 = inM[7]; // M(2,1)
-
- const double m02 = inM[2]; // M(0,2)
- const double m12 = inM[5]; // M(1,2)
- const double m22 = inM[8]; // M(2,2)
-
- // ----------------------------------------------------------------------
- // 1) Compute determinant:
- // det(M) = m00*(m11*m22 - m12*m21)
- // - m01*(m10*m22 - m12*m20)
- // + m02*(m10*m21 - m11*m20)
- // ----------------------------------------------------------------------
- double det = m00*(m11*m22 - m12*m21)
- - m01*(m10*m22 - m12*m20)
- + m02*(m10*m21 - m11*m20);
-
- if (std::fabs(det) < 1e-14) {
- throw std::runtime_error(
- "inverse_3x3_cofactor_custom: matrix is near-singular."
- );
- }
- double invDet = 1.0 / det;
-
- // ----------------------------------------------------------------------
- // 2) Cofactor / adjugate formula:
- //
- // inv(0,0) = (m11*m22 - m12*m21)
- // inv(0,1) = -(m01*m22 - m02*m21)
- // inv(0,2) = (m01*m12 - m02*m11)
- //
- // inv(1,0) = -(m10*m22 - m12*m20)
- // inv(1,1) = (m00*m22 - m02*m20)
- // inv(1,2) = -(m00*m12 - m02*m10)
- //
- // inv(2,0) = (m10*m21 - m11*m20)
- // inv(2,1) = -(m00*m21 - m01*m20)
- // inv(2,2) = (m00*m11 - m01*m10)
- //
- // Multiply each cofactor by invDet.
- // ----------------------------------------------------------------------
- double i00 = (m11*m22 - m12*m21) * invDet;
- double i01 = -(m01*m22 - m02*m21) * invDet;
- double i02 = (m01*m12 - m02*m11) * invDet;
-
- double i10 = -(m10*m22 - m12*m20) * invDet;
- double i11 = (m00*m22 - m02*m20) * invDet;
- double i12 = -(m00*m12 - m02*m10) * invDet;
-
- double i20 = (m10*m21 - m11*m20) * invDet;
- double i21 = -(m00*m21 - m01*m20) * invDet;
- double i22 = (m00*m11 - m01*m10) * invDet;
-
- // ----------------------------------------------------------------------
- // 3) Store the inverse back into outM in the same layout:
- // inv(r,c) => outM[r*3 + c]
- // ----------------------------------------------------------------------
- outM[0] = i00; // inv(0,0)
- outM[1] = i01; // inv(0,1)
- outM[2] = i02; // inv(0,2)
-
- outM[3] = i10; // inv(1,0)
- outM[4] = i11; // inv(1,1)
- outM[5] = i12; // inv(1,2)
-
- outM[6] = i20; // inv(2,0)
- outM[7] = i21; // inv(2,1)
- outM[8] = i22; // inv(2,2)
-}
-
- /**
- * Naive approach to build neighbor lists.
- */
- void calc_naive(Structure &st);
-
- /**
- * Binning-based approach. If any dimension < cutoff, fallback to naive.
- */
- void calc_binned(Structure &st);
-
- /**
- * For naive approach, compute ±N to consider for image shifts.
- */
- void num_shifts(Structure &st, int N[3]);
-
-public:
- /**
- * Construct with cutoff from config("RCUTMAX").
- */
- NNFinder(Config &config);
-
- /**
- * Build nearest neighbors for all atoms in one Structure.
- * Uses binned approach if possible, else naive fallback.
- */
- void calc(Structure &st);
-
- /**
- * Build nearest neighbors for each structure in a DB.
- */
- void calc(StructureDB &stdb);
-};
-}
-}
-
-#endif // NN_FINDER_H
-
-
diff --git a/include/tadah/mlip/normaliser.h b/include/tadah/mlip/normaliser.h
index f78a1c7a0084df1ae13defdce925042a623f8143..1e0e3bb76cc96a74ab0642edc8f07a903b51296a 100644
--- a/include/tadah/mlip/normaliser.h
+++ b/include/tadah/mlip/normaliser.h
@@ -11,18 +11,18 @@
namespace tadah {
namespace mlip {
-class Normaliser: public Normaliser_Core {
+class Normaliser: public tadah::core::Normaliser_Core {
public:
- using Normaliser_Core::normalise;
+ using tadah::core::Normaliser_Core::normalise;
Normaliser ():
- Normaliser_Core()
+ tadah::core::Normaliser_Core()
{};
- Normaliser (Config &c):
- Normaliser_Core(c)
+ Normaliser (tadah::core::Config &c):
+ tadah::core::Normaliser_Core(c)
{};
- Normaliser(Config &c,StDescriptorsDB &st_desc_db):
- Normaliser_Core(c)
+ Normaliser(tadah::core::Config &c,StDescriptorsDB &st_desc_db):
+ tadah::core::Normaliser_Core(c)
{
learn(st_desc_db);
@@ -38,13 +38,13 @@ class Normaliser: public Normaliser_Core {
// prep containers
size_t dim = st_desc_db(0).dim();
- mean = v_type(dim);
- std_dev = v_type(dim);
+ mean = tadah::core::v_type(dim);
+ std_dev = tadah::core::v_type(dim);
//std::cout << "COLS NORM: " << nn << std::endl;
// compute mean and st_dev
for (size_t d=0; d<dim; ++d) {
- t_type v(n);
+ tadah::core::t_type v(n);
size_t b=0;
for (size_t s=0; s<st_desc_db.size(); ++s) {
//for (const auto &aed:st_desc_db(s).aed) {
@@ -73,6 +73,7 @@ class Normaliser: public Normaliser_Core {
}
}
+ using tadah::core::operator<<;
if (verbose) std::cout << std::endl << "NORMALISER STDEV : " << std_dev << std::endl;
if (verbose) std::cout << std::endl << "NORMALISER MEAN :" << mean << std::endl;
diff --git a/include/tadah/mlip/output/output.h b/include/tadah/mlip/output/output.h
index 1a0885035b3772d629927af266a9e55a460917b1..690d0bd0012a0c09feccb87c6b9bd7a26e0a20a3 100644
--- a/include/tadah/mlip/output/output.h
+++ b/include/tadah/mlip/output/output.h
@@ -13,7 +13,7 @@
class Output {
private:
- Config &pot_config;
+ tadah::core::Config &pot_config;
bool pred_err_bool;
//size_t sigfig;
//size_t nw;
@@ -22,11 +22,11 @@ class Output {
size_t emptyspace=5;
size_t nw=emptyspace+fmtw; // width of numeric column
public:
- Output (Config &c, bool peb):
+ Output (tadah::core::Config &c, bool peb):
pot_config(c),
pred_err_bool(peb)
{}
- void print_train_unc(t_type &weights, t_type &unc) {
+ void print_train_unc(tadah::core::t_type &weights, tadah::core::t_type &unc) {
std::string index_str="Index";
size_t iwe=(size_t)unc.size() > max_number(index_str.size()) ?
@@ -45,7 +45,7 @@ class Output {
out_unc.close();
}
- void print_predict_all(StructureDB &stdb, StructureDB &stpred, aed_type & predicted_error) {
+ void print_predict_all(StructureDB &stdb, StructureDB &stpred, tadah::core::aed_type & predicted_error) {
std::ofstream out_error("error.pred");
std::ofstream out_energy("energy.pred");
std::ofstream out_force("forces.pred");
@@ -143,8 +143,8 @@ class Output {
if (newline)
out_force << std::endl << std::endl;
for (size_t a=0; a<stdb(i).natoms(); ++a) {
- Vec3d &f = stdb(i).atoms[a].force;
- Vec3d &fp = stpred(i).atoms[a].force;
+ tadah::core::Vec3d &f = stdb(i).atoms[a].force;
+ tadah::core::Vec3d &fp = stpred(i).atoms[a].force;
for (size_t k=0; k<3; ++k) {
out_force << std::setw(iwf) << fidx++ << std::setw(nw)
<< f[k] << std::setw(nw) << fp[k];
@@ -163,8 +163,8 @@ class Output {
double max_serr=0.0;
if (newline)
out_stress << std::endl << std::endl;
- stress_type &s = stdb(i).stress;
- stress_type &sp = stpred(i).stress;
+ tadah::core::stress_type &s = stdb(i).stress;
+ tadah::core::stress_type &sp = stpred(i).stress;
for (size_t x=0; x<3; ++x) {
for (size_t y=x; y<3; ++y) {
out_stress << std::setw(iws) << sidx++ << std::setw(nw) << s(x,y)
diff --git a/include/tadah/mlip/st_descriptors.h b/include/tadah/mlip/st_descriptors.h
index a3206efdd3300dcc5f37df1186cb0ac627eddfb8..c758449a4074cad34fd8af27d7b5749d936284ab 100644
--- a/include/tadah/mlip/st_descriptors.h
+++ b/include/tadah/mlip/st_descriptors.h
@@ -20,11 +20,11 @@ namespace mlip {
* To fully initialise this object it needs to know:
* - number of atoms in a structure (from Structure)
* - number of nn for every atom (from Structure)
- * - dimension of the descriptor vector (Config)
- * - whether force and stress is being calculated (from Config)
+ * - dimension of the descriptor vector (tadah::core::Config)
+ * - whether force and stress is being calculated (from tadah::core::Config)
*
* \note
- * Required Config keys:
+ * Required tadah::core::Config keys:
* \ref FORCE \ref STRESS.
* \ref INTERNAL_KEY \ref DSIZE.
*
@@ -35,9 +35,9 @@ struct StDescriptors {
*
* Requires:
* - Structure st to have NN calculated
- * - Config c to contain keys: \ref DSIZE, \ref FORCE, \ref STRESS
+ * - tadah::core::Config c to contain keys: \ref DSIZE, \ref FORCE, \ref STRESS
*/
- StDescriptors(const Structure &s, const Config &c);
+ StDescriptors(const Structure &s, const tadah::core::Config &c);
/** Default constructor. Object is left uninitialised */
StDescriptors();
@@ -50,14 +50,14 @@ struct StDescriptors {
* Dimensions:
* [number of atoms, descriptor size]
*/
- aeds_type2 aeds;
+ tadah::core::aeds_type2 aeds;
/** FD for all atoms
*
* Dimensions:
* [number of atoms, number of atom NN, descriptor size, 3]
*/
- std::vector<std::vector<fd_type>> fd;
+ std::vector<std::vector<tadah::core::fd_type>> fd;
/** SD for a Structure
*
@@ -69,11 +69,11 @@ struct StDescriptors {
/** True if descriptors are normalised */
bool normalised=false;
- rhos_type rhos;
+ tadah::core::rhos_type rhos;
- aed_type & get_aed(const size_t i);
- const aed_type &get_aed(const size_t i) const;
- rho_type& get_rho(const size_t i);
+ tadah::core::aed_type & get_aed(const size_t i);
+ const tadah::core::aed_type &get_aed(const size_t i) const;
+ tadah::core::rho_type& get_rho(const size_t i);
size_t naed() const;
size_t dim() const;
diff --git a/include/tadah/mlip/st_descriptors_db.h b/include/tadah/mlip/st_descriptors_db.h
index d9aa5ed612d084fb91f1e20ceb4b19ac43583677..005d95db0ef1069d1a8676a8c04eccb071012881 100644
--- a/include/tadah/mlip/st_descriptors_db.h
+++ b/include/tadah/mlip/st_descriptors_db.h
@@ -11,7 +11,7 @@ namespace mlip {
/** \brief Container for StDescriptors.
*
* \note
- * Required Config keys:
+ * Required tadah::core::Config keys:
* \ref FORCE \ref STRESS
* \ref INTERNAL_KEY \ref DSIZE.
*
@@ -24,10 +24,10 @@ struct StDescriptorsDB {
*
* Requires:
* - StructureDB st to have all nearest neighbours calculated
- * - Config to contain keys: \ref FORCE, \ref STRESS
+ * - tadah::core::Config to contain keys: \ref FORCE, \ref STRESS
* and \ref INTERNAL_KEY \ref DSIZE
*/
- StDescriptorsDB(const StructureDB &stdb, Config &config);
+ StDescriptorsDB(const StructureDB &stdb, tadah::core::Config &config);
/** Return reference to the n-th StDescriptors
* stored by this object
diff --git a/include/tadah/mlip/structure.h b/include/tadah/mlip/structure.h
index 908f8180179ce5fc8822df5b85754d3e4321eb0c..89063e039e9814e8e264aa6b02abd89345d39b5b 100644
--- a/include/tadah/mlip/structure.h
+++ b/include/tadah/mlip/structure.h
@@ -8,8 +8,9 @@
#include <stdexcept>
#include <iterator> // for std::forward_iterator_tag
#include <set>
+// Do not include atom.h
-// HPC types (Matrix3d, stress_type) come from here:
+// HPC types (tadah::core::Matrix3d, tadah::core::stress_type) come from here:
#include <tadah/core/periodic_table.h>
#include <tadah/core/core_types.h>
#include <tadah/core/vec3d_soa_view.h>
@@ -38,8 +39,8 @@ public:
std::string label;
double energy = 0.0;
- Matrix3d cell;
- stress_type stress;
+ tadah::core::Matrix3d cell;
+ tadah::core::stress_type stress;
double T = 0.0;
double eweight = 1.0, fweight = 1.0, sweight = 1.0;
@@ -85,11 +86,11 @@ public:
const int& atomicNumber(size_t i) const;
// SoA views for HPC (non-const)
- Vec3dSoAView positionView(size_t i);
- Vec3dSoAView forceView(size_t i);
+ tadah::core::Vec3dSoAView positionView(size_t i);
+ tadah::core::Vec3dSoAView forceView(size_t i);
- Vec3dSoAConstView positionView(size_t i) const;
- Vec3dSoAConstView forceView(size_t i) const;
+ tadah::core::Vec3dSoAConstView positionView(size_t i) const;
+ tadah::core::Vec3dSoAConstView forceView(size_t i) const;
/**
* @brief Converts a global index to a local index within this structure.
diff --git a/include/tadah/mlip/structure_db.h b/include/tadah/mlip/structure_db.h
index 89090bacd9f70b37326e3d353baea3734c68f672..3900b5d621ae14101c0672d7d1e475dc1607df73 100644
--- a/include/tadah/mlip/structure_db.h
+++ b/include/tadah/mlip/structure_db.h
@@ -2,163 +2,608 @@
#ifndef TADAH_MLIP_STRUCTURE_DB_H
#define TADAH_MLIP_STRUCTURE_DB_H
-#include <vector>
-#include <string>
+#include <cctype>
#include <fstream>
#include <iostream>
-#include <sstream>
#include <set>
-#include <cctype>
+#include <sstream>
#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <vector>
#ifdef _OPENMP
-#include <omp.h> // Optional if parallelizing
+#include <omp.h> // Optional if parallelizing
#endif
-#include <tadah/core/config.h> // tadah::core::Config
-#include <tadah/mlip/structure.h> // tadah::mlip::Structure
+#include <tadah/core/config.h> // tadah::core::Config
+#include <tadah/mlip/structure.h> // tadah::mlip::Structure
namespace tadah {
namespace mlip {
/**
- * @brief Holds a set of structures loaded from multiple files, using a two-pass approach.
+ * @brief Holds a set of structures loaded from multiple files, using a two-pass
+ * approach.
*
* The approach:
- * 1) Parse how many structures and how many atoms each structure has (pass1).
- * 2) Allocate SoA arrays for HPC storage (X_,Y_,Z_, etc.).
- * 3) Re-read each file to fill data for each structure (pass2).
+ * 1) Pass1: Parse how many structures and how many atoms each structure has.
+ * 2) Allocate SoA arrays for HPC storage (X_, Y_, Z_, FX_, FY_, FZ_,
+ * atomicNumber_). 3) Pass2: Re-read each file to fill data for each structure.
*
- * Each file can contain multiple structures separated by blank lines.
- * Each structure has a 9-line header plus N atom lines. The HPC arrays store
- * positions, forces, and atomic numbers in contiguous format.
+ * Each file can contain multiple structures separated by blank lines.
+ * Each structure has a 9-line header plus N lines describing atom data.
+ * The HPC arrays store positions, forces, and atomic numbers in contiguous
+ * format.
*/
class StructureDB {
public:
-
- /**
- * @brief Stores pass1 results such as how many structures are in each file,
- * how many atoms each structure contains, and which files are included.
- *
- * Permits loading a certain subset of structures from large input files by
- * specifying which ranges are included.
- */
+ /**
+ * @brief Pass1Info stores results of the first pass: how many structures, how
+ * many atoms per structure, etc.
+ */
struct Pass1Info {
Pass1Info() = default;
- Pass1Info(const std::vector<std::string> &files,
- const std::vector<std::size_t> &nStructuresPerFile,
- const std::vector<std::vector<std::size_t>> &atomsPerStructurePerFile)
- : files(files), nStructuresPerFile(nStructuresPerFile), atomsPerStructurePerFile(atomsPerStructurePerFile)
- {
+ Pass1Info(
+ const std::vector<std::string> &files,
+ const std::vector<std::size_t> &nStructuresPerFile,
+ const std::vector<std::vector<std::size_t>> &atomsPerStructurePerFile)
+ : files(files), nStructuresPerFile(nStructuresPerFile),
+ atomsPerStructurePerFile(atomsPerStructurePerFile) {
firstStructureIndexPerFile.resize(files.size(), 0);
}
- Pass1Info(const std::vector<std::string> &files,
- const std::vector<std::size_t> &nStructuresPerFile,
- const std::vector<std::vector<std::size_t>> &atomsPerStructurePerFile,
- const std::vector<std::size_t> &firstStructureIndexPerFile):
- files(files), nStructuresPerFile(nStructuresPerFile),
- atomsPerStructurePerFile(atomsPerStructurePerFile),
- firstStructureIndexPerFile(firstStructureIndexPerFile) {}
-
+ Pass1Info(
+ const std::vector<std::string> &files,
+ const std::vector<std::size_t> &nStructuresPerFile,
+ const std::vector<std::vector<std::size_t>> &atomsPerStructurePerFile,
+ const std::vector<std::size_t> &firstStructureIndexPerFile)
+ : files(files), nStructuresPerFile(nStructuresPerFile),
+ atomsPerStructurePerFile(atomsPerStructurePerFile),
+ firstStructureIndexPerFile(firstStructureIndexPerFile) {}
+
std::vector<std::string> files;
std::vector<std::size_t> nStructuresPerFile;
std::vector<std::size_t> firstStructureIndexPerFile;
std::vector<std::vector<std::size_t>> atomsPerStructurePerFile;
};
-
+ /**
+ * @brief Returns a set of unique atomic numbers across all structures.
+ */
std::set<int> get_unique_elements() const {
std::set<int> s;
- for (const auto & atomic_num : atomicNumber_) {
+ for (const auto &atomic_num : atomicNumber_) {
s.insert(atomic_num);
}
return s;
}
+private:
+ // -------------------------------------------------------------------------
+ // Constants controlling file parsing
+ // -------------------------------------------------------------------------
+ static constexpr std::size_t HEADERSIZE_ =
+ 9; ///< Each structure block has 9 header lines
+ static constexpr std::size_t BUFSIZE_ =
+ (100ULL << 20); ///< For potential large-file buffering
+ static constexpr int TANTALUM_Z_ =
+ 73; ///< Example usage if "Ta" is encountered
+
+public:
+ // -------------------------------------------------------------------------
+ // HPC SoA storage
+ // -------------------------------------------------------------------------
+ std::vector<double> X_, Y_, Z_; ///< Positions
+ std::vector<double> FX_, FY_, FZ_; ///< Forces
+ std::vector<int> atomicNumber_;
private:
- /**
- * @brief HEADERSIZE_ references the 9 lines forming a header in each structure file block.
- */
- static constexpr std::size_t HEADERSIZE_ = 9;
+ // -------------------------------------------------------------------------
+ // Internal container of structures
+ // -------------------------------------------------------------------------
+ std::vector<Structure> structures_;
+ std::size_t usedAtoms_ = 0;
- /**
- * @brief BUFSIZE_ sets a 100MB read buffer to accelerate file I/O for large data.
- */
- static constexpr std::size_t BUFSIZE_ = (100ULL << 20);
+ /**
+ * @brief structureAtomCountsPerFile_[f] stores the number of atoms in each
+ * structure of file f.
+ */
+ std::vector<std::vector<std::size_t>> structureAtomCountsPerFile_;
- /**
- * @brief TANTALUM_Z_ used as a dummy atomic number for "Ta" examples.
- */
- static constexpr int TANTALUM_Z_ = 73;
+ /**
+ * @brief fileStructureOffset_[f] indicates the first structure index in
+ * "structures_" that belongs to file f.
+ */
+ std::vector<std::size_t> fileStructureOffset_;
public:
- std::vector<double> X_, Y_, Z_;
- std::vector<double> FX_, FY_, FZ_;
- std::vector<int> atomicNumber_;
+ /**
+ * @brief Constructs from a config object: expects config entries listing the
+ * input files.
+ *
+ * Implementation:
+ * - Reads the list of input files from config.
+ * - For each file, calls pass1 to determine how many structures and atoms.
+ * - Allocates SoA arrays.
+ * - Invokes pass2 to fill data for each structure.
+ *
+ * Example usage:
+ * @code
+ * tadah::core::Config config("myconfig.json");
+ * StructureDB db(config);
+ * std::cout << db.summary() << std::endl;
+ * @endcode
+ */
+ inline explicit StructureDB(tadah::core::Config &config) {
+ // Reason: The config approach is application-specific.
+ // This code attempts to retrieve a vector of input filenames
+ // from a key "files", then runs the two-pass procedure.
+
+ std::vector<std::string> files;
+ try {
+ config.get<std::vector<std::string>>("DBFILE", files);
+ } catch (...) {
+ throw std::runtime_error(
+ "StructureDB: 'files' not found in config or invalid type");
+ }
+
+ // pass1
+ std::vector<std::size_t> nStructuresPerFile;
+ std::vector<std::vector<std::size_t>> atomsPerStructurePerFile;
+ nStructuresPerFile.reserve(files.size());
+ atomsPerStructurePerFile.reserve(files.size());
+
+ for (const auto &f : files) {
+ // gather per-file structure-atom counts
+ auto structureAtomCounts = pass1ParseStructures(f);
+ nStructuresPerFile.push_back(structureAtomCounts.size());
+ atomsPerStructurePerFile.push_back(structureAtomCounts);
+ }
+
+ // Build pass1 info
+ Pass1Info pinfo(files, nStructuresPerFile, atomsPerStructurePerFile);
+ storedPass1_ = pinfo; // Store the pass1 info
+ // Allocate HPC arrays
+ // Reason: Sum up total atoms across all structures, then resize vectors.
+ std::size_t totalStructures = 0;
+ std::size_t totalAtoms = 0;
+ for (std::size_t i = 0; i < files.size(); i++) {
+ totalStructures += pinfo.nStructuresPerFile[i];
+ for (auto count : pinfo.atomsPerStructurePerFile[i]) {
+ totalAtoms += count;
+ }
+ }
+ X_.resize(totalAtoms);
+ Y_.resize(totalAtoms);
+ Z_.resize(totalAtoms);
+ FX_.resize(totalAtoms);
+ FY_.resize(totalAtoms);
+ FZ_.resize(totalAtoms);
+ atomicNumber_.resize(totalAtoms);
+
+ // Prepare structure offsets per file
+ // Reason: track how many structures came before file f
+ fileStructureOffset_.resize(files.size(), 0);
+ {
+ std::size_t offsetSoFar = 0;
+ for (std::size_t fIndex = 0; fIndex < files.size(); fIndex++) {
+ fileStructureOffset_[fIndex] = offsetSoFar;
+ offsetSoFar += pinfo.nStructuresPerFile[fIndex];
+ }
+ }
+
+ // Create actual Structure objects in structures_
+ structures_.reserve(totalStructures);
+ usedAtoms_ = 0;
+ for (std::size_t fIndex = 0; fIndex < files.size(); fIndex++) {
+ // For each structure in file fIndex
+ for (std::size_t sIndex = 0; sIndex < pinfo.nStructuresPerFile[fIndex];
+ sIndex++) {
+ // This structure will begin at usedAtoms_ in HPC arrays
+ auto countAtoms = pinfo.atomsPerStructurePerFile[fIndex][sIndex];
+ structures_.emplace_back(this, usedAtoms_, countAtoms);
+ usedAtoms_ += countAtoms;
+ }
+ }
+
+ // pass2 fill
+ for (std::size_t fIndex = 0; fIndex < files.size(); fIndex++) {
+ pass2FillStructures(fIndex, fileStructureOffset_[fIndex], files[fIndex]);
+ }
+ }
+
+ /**
+ * @brief Constructs from external pass1 info. Expects the caller to have
+ * precomputed how many structures and atoms come from each file.
+ *
+ * Example usage:
+ * @code
+ * // Suppose pass1 was done externally
+ * StructureDB::Pass1Info info(...);
+ * StructureDB db(info);
+ * @endcode
+ */
+ inline StructureDB(const Pass1Info &pinfo) {
+ storedPass1_ = pinfo; // Store the pass1 info
+ // Summation of total structures and atoms
+ std::size_t totalStructures = 0;
+ std::size_t totalAtoms = 0;
+ for (std::size_t i = 0; i < pinfo.files.size(); i++) {
+ totalStructures += pinfo.nStructuresPerFile[i];
+ for (auto count : pinfo.atomsPerStructurePerFile[i]) {
+ totalAtoms += count;
+ }
+ }
+
+ // Allocate SoA arrays
+ X_.resize(totalAtoms);
+ Y_.resize(totalAtoms);
+ Z_.resize(totalAtoms);
+ FX_.resize(totalAtoms);
+ FY_.resize(totalAtoms);
+ FZ_.resize(totalAtoms);
+ atomicNumber_.resize(totalAtoms);
+
+ // Prepare structure offsets per file
+ fileStructureOffset_.resize(pinfo.files.size(), 0);
+ {
+ std::size_t offsetSoFar = 0;
+ for (std::size_t fIndex = 0; fIndex < pinfo.files.size(); fIndex++) {
+ fileStructureOffset_[fIndex] = offsetSoFar;
+ offsetSoFar += pinfo.nStructuresPerFile[fIndex];
+ }
+ }
+
+ // Create the structure objects
+ structures_.reserve(totalStructures);
+ usedAtoms_ = 0;
+ for (std::size_t fIndex = 0; fIndex < pinfo.files.size(); fIndex++) {
+ for (std::size_t sIndex = 0; sIndex < pinfo.nStructuresPerFile[fIndex];
+ sIndex++) {
+ auto countAtoms = pinfo.atomsPerStructurePerFile[fIndex][sIndex];
+ structures_.emplace_back(this, usedAtoms_, countAtoms);
+ usedAtoms_ += countAtoms;
+ }
+ }
+
+ // pass2 approach: re-read each file for actual data
+ // The user might choose to do it manually or replicate pass2 logic
+ // externally.
+ for (std::size_t fIndex = 0; fIndex < pinfo.files.size(); fIndex++) {
+ pass2FillStructures(fIndex, fileStructureOffset_[fIndex],
+ pinfo.files[fIndex]);
+ }
+ }
+
+ /**
+ * @brief Constructs from external pass1 info, optionally skipping pass2.
+ *
+ * The pass1 info determines the number of total structures and total atoms.
+ * HPC arrays are allocated and the Structure instances are created but not
+ * necessarily filled from file. This makes a DB with the "shape" of pinfo
+ * but without data if skipPass2 is true.
+ *
+ * Example usage:
+ * @code
+ * // Assume some existing DB with pass1 info
+ * StructureDB db(...); // does pass2 fill
+ * auto pinfo = db.exportPass1Info();
+ *
+ * // Build an "empty" copy with same structure sizes but no pass2
+ * StructureDB db_empty(pinfo, true);
+ * std::cout << db_empty.summary() << "\n"; // same # of structures
+ * @endcode
+ *
+ * @param pinfo the pass1 info describing how many files, structures, atoms
+ * @param skipPass2 if true, pass2FillStructures is skipped.
+ */
+ inline StructureDB(const Pass1Info &pinfo, bool skipPass2) {
+ // Store pass1 info for future reference
+ storedPass1_ = pinfo;
+
+ // Summation of total structures and atoms
+ std::size_t totalStructures = 0;
+ std::size_t totalAtoms = 0;
+ for (std::size_t i = 0; i < pinfo.files.size(); i++) {
+ totalStructures += pinfo.nStructuresPerFile[i];
+ for (auto count : pinfo.atomsPerStructurePerFile[i]) {
+ totalAtoms += count;
+ }
+ }
+
+ // Allocate SoA arrays
+ X_.resize(totalAtoms);
+ Y_.resize(totalAtoms);
+ Z_.resize(totalAtoms);
+ FX_.resize(totalAtoms);
+ FY_.resize(totalAtoms);
+ FZ_.resize(totalAtoms);
+ atomicNumber_.resize(totalAtoms);
+
+ // Prepare structure offsets per file
+ fileStructureOffset_.resize(pinfo.files.size(), 0);
+ {
+ std::size_t offsetSoFar = 0;
+ for (std::size_t fIndex = 0; fIndex < pinfo.files.size(); fIndex++) {
+ fileStructureOffset_[fIndex] = offsetSoFar;
+ offsetSoFar += pinfo.nStructuresPerFile[fIndex];
+ }
+ }
+
+ // Create the structure objects in structures_
+ structures_.reserve(totalStructures);
+ usedAtoms_ = 0;
+ for (std::size_t fIndex = 0; fIndex < pinfo.files.size(); fIndex++) {
+ // Store per-file atom counts for the top-level queries
+ structureAtomCountsPerFile_.push_back(
+ pinfo.atomsPerStructurePerFile[fIndex]);
+
+ // Create each structure
+ for (std::size_t sIndex = 0; sIndex < pinfo.nStructuresPerFile[fIndex];
+ sIndex++) {
+ auto countAtoms = pinfo.atomsPerStructurePerFile[fIndex][sIndex];
+ structures_.emplace_back(this, usedAtoms_, countAtoms);
+ usedAtoms_ += countAtoms;
+ }
+ }
+
+ // Perform pass2 only if skipPass2 == false
+ // If skipPass2 == true, HPC arrays remain uninitialized
+ if (!skipPass2) {
+ for (std::size_t fIndex = 0; fIndex < pinfo.files.size(); fIndex++) {
+ pass2FillStructures(fIndex, fileStructureOffset_[fIndex],
+ pinfo.files[fIndex]);
+ }
+ }
+ }
+
+ /**
+ * @brief Returns total number of structures across all files.
+ */
+ inline std::size_t size() const { return structures_.size(); }
+
+ /**
+ * @brief Returns how many files were processed.
+ */
+ inline std::size_t nFiles() const {
+ return structureAtomCountsPerFile_.size();
+ // Note: This is non-empty after pass1 or the second constructor's
+ // population
+ }
+
+ /**
+ * @brief Returns how many structures were found in file f.
+ * @param f index of the file
+ */
+ inline std::size_t nStructuresInFile(std::size_t f) const {
+ if (f >= structureAtomCountsPerFile_.size()) {
+ throw std::out_of_range(
+ "StructureDB::nStructuresInFile: file index out of range");
+ }
+ return structureAtomCountsPerFile_[f].size();
+ }
+
+ /**
+ * @brief Access to the i-th structure overall.
+ */
+ inline Structure &operator()(std::size_t i) {
+ if (i >= structures_.size()) {
+ throw std::out_of_range("StructureDB::operator()(i): index out of range");
+ }
+ return structures_[i];
+ }
+ inline const Structure &operator()(std::size_t i) const {
+ if (i >= structures_.size()) {
+ throw std::out_of_range(
+ "StructureDB::operator()(i) const: index out of range");
+ }
+ return structures_[i];
+ }
+
+ /**
+ * @brief Access to structure j in file f.
+ */
+ inline Structure &operator()(std::size_t f, std::size_t j) {
+ if (f >= fileStructureOffset_.size()) {
+ throw std::out_of_range(
+ "StructureDB::operator()(f, j): file index out of range");
+ }
+ std::size_t i = fileStructureOffset_[f] + j;
+ return (*this)(i);
+ }
+ inline const Structure &operator()(std::size_t f, std::size_t j) const {
+ if (f >= fileStructureOffset_.size()) {
+ throw std::out_of_range(
+ "StructureDB::operator()(f, j) const: file index out of range");
+ }
+ std::size_t i = fileStructureOffset_[f] + j;
+ return (*this)(i);
+ }
+
+ /**
+ * @brief Summary of how many structures are loaded in total.
+ */
+ inline std::string summary() const {
+ std::ostringstream oss;
+ oss << "StructureDB: " << size() << " structures total\n";
+ return oss.str();
+ }
private:
- std::vector<Structure> structures_;
- std::size_t usedAtoms_ = 0;
+ /**
+ * @brief Pass1: parses each file to identify how many structures
+ * and how many atoms each structure contains.
+ */
+ inline std::vector<std::size_t>
+ pass1ParseStructures(const std::string &filename) const {
+ // Reason: pass1 only inspects the file line-by-line, counting
+ // how many lines belong to each structure until blank lines appear.
+ // Each structure block has 9 lines of header + N lines of atoms.
+ // The total lines from the block end at the next blank line.
- std::vector<std::vector<std::size_t>> structureAtomCountsPerFile_;
- std::vector<std::size_t> fileStructureOffset_;
+ std::ifstream ifs(filename);
+ if (!ifs.is_open()) {
+ throw std::runtime_error(
+ "StructureDB::pass1ParseStructures: Cannot open file " + filename);
+ }
-public:
- /**
- * @brief Constructs by reading config for DB files, performs two-pass fill (if desired).
- */
- inline explicit StructureDB(Config &config);
+ std::vector<std::size_t> atomsPerStructure;
+ while (true) {
+ // Attempt to read 9 header lines for a new structure
+ std::string line;
+ // Skip potential blanks
+ while (std::getline(ifs, line)) {
+ if (!isBlankLine(line)) {
+ // Puts the line back into the stream
+ ifs.seekg(-static_cast<std::streamoff>(line.size() + 1),
+ std::ios::cur);
+ break;
+ }
+ }
+ if (!ifs.good())
+ break;
- /**
- * @brief Constructs from external pass1 info, then allocates HPC arrays accordingly.
- */
- inline StructureDB(const Pass1Info &pinfo);
+ // Read and discard 9 lines for the structure header if possible
+ bool headerOk = true;
+ for (std::size_t i = 0; i < HEADERSIZE_; i++) {
+ if (!std::getline(ifs, line)) {
+ headerOk = false;
+ break;
+ }
+ }
+ if (!headerOk)
+ break; // Reached EOF inside header
- // Returns total number of structures across all files
- inline std::size_t size() const { return structures_.size(); }
+ // Now read atom lines until a blank line or EOF
+ std::size_t atomCount = 0;
+ std::streampos pos;
+ while (true) {
+ pos = ifs.tellg();
+ if (!std::getline(ifs, line)) {
+ break; // EOF
+ }
+ if (isBlankLine(line)) {
+ // Structure block ended
+ break;
+ }
+ // Attempt to parse an atom line
+ std::stringstream tmp(line);
+ int zVal;
+ double xx, yy, zz, fxx, fyy, fzz;
+ tmp >> zVal >> xx >> yy >> zz >> fxx >> fyy >> fzz;
+ if (tmp.fail()) {
+ // Not a valid atom line => roll back, break
+ ifs.seekg(pos);
+ break;
+ }
+ atomCount++;
+ }
+ atomsPerStructure.push_back(atomCount);
- // Returns how many files were processed
- inline std::size_t nFiles() const { return structureAtomCountsPerFile_.size(); }
+ // That structure is done, next iteration might read next structure or
+ // end.
+ if (!ifs.good() && ifs.eof()) {
+ break;
+ }
+ }
+ // Close and return
+ ifs.close();
+ return atomsPerStructure;
+ }
- // Returns how many structures are found in file f
- inline std::size_t nStructuresInFile(std::size_t f) const;
+ /**
+ * @brief Pass2: re-opens each file, re-reads each structure’s header & atom
+ * lines, and fills the HPC SoA arrays for every structure in that file.
+ *
+ * @param fIndex index of the file
+ * @param offsetS starting structure index in "structures_" for this file
+ * @param ifile the name of the file to parse
+ */
+ inline void pass2FillStructures(std::size_t fIndex, std::size_t offsetS,
+ const std::string &ifile) {
+ std::ifstream ifs(ifile);
+ if (!ifs.is_open()) {
+ throw std::runtime_error(
+ "StructureDB::pass2FillStructures: Cannot open file " + ifile);
+ }
+ // Reason: A read buffer can be set if desired for large files (not
+ // mandatory).
+ ifs.rdbuf()->pubsetbuf(nullptr, BUFSIZE_);
- // Access to structures by overall index or by (file, structureIndexInFile)
- inline Structure& operator()(std::size_t i);
- inline const Structure& operator()(std::size_t i) const;
+ std::size_t sCount = 0; // how many structures from this file so far
+ while (true) {
+ // Skip blank lines
+ std::string line;
+ while (std::getline(ifs, line)) {
+ if (!isBlankLine(line)) {
+ ifs.seekg(-static_cast<std::streamoff>(line.size() + 1),
+ std::ios::cur);
+ break;
+ }
+ }
+ if (!ifs.good())
+ break;
+ if (sCount >= structureAtomCountsPerFile_[fIndex].size())
+ break;
- inline Structure& operator()(std::size_t f, std::size_t j);
- inline const Structure& operator()(std::size_t f, std::size_t j) const;
+ // The next structure to fill is structures_[offsetS + sCount]
+ // Use Structure::read() method to parse the block
+ auto &struc = structures_[offsetS + sCount];
+ struc.read(ifs);
+ sCount++;
+ if (!ifs.good() && ifs.eof()) {
+ break;
+ }
+ }
+ ifs.close();
+ }
- // Summary for debug
- inline std::string summary() const {
- std::ostringstream oss;
- oss << "StructureDB: " << size() << " structures total\n";
- return oss.str();
+ /**
+ * @brief Helper to identify blank or whitespace-only lines.
+ */
+ inline bool isBlankLine(const std::string &line) const {
+ // Checks if all characters are space or tab
+ for (char c : line) {
+ if (!std::isspace(static_cast<unsigned char>(c))) {
+ return false;
+ }
}
+ return true;
+ }
+
+ /**
+ * @brief Converts a chemical symbol to an atomic number.
+ * Returns TANTALUM_Z_ for "Ta", else zeros out by default.
+ * (Example usage only; the real code may implement a more complete
+ * logic.)
+ */
+ inline int toAtomicNumber(const std::string &sym) const {
+ // Reason: Minimal placeholder for demonstration.
+ // Real code might reference a more complete periodic table.
+ if (sym == "Ta")
+ return TANTALUM_Z_;
+ return 0;
+ }
private:
- /**
- * @brief Pass1 parses each file quickly to identify how many structures
- * and how many lines each structure contains.
- */
- inline std::vector<std::size_t> pass1ParseStructures(const std::string &filename) const;
-
- /**
- * @brief Pass2 re-opens each file, skips the 9-line header, and reads
- * the atom lines, storing them in HPC SoA arrays.
- */
- inline void pass2FillStructures(std::size_t fIndex,
- std::size_t offsetS,
- const std::string &ifile);
-
- inline bool isBlankLine(const std::string& line) const;
-
- /**
- * @brief Converts a symbol to an atomic number. Returns TANTALUM_Z_ if "Ta," else returns 0.
- */
- inline int toAtomicNumber(const std::string &sym) const;
+ /**
+ * @brief Stores the pass1 info used by this DB, allowing easy export later.
+ */
+ Pass1Info storedPass1_;
+
+public:
+ /**
+ * @brief Exports the internally stored Pass1Info to replicate the structure
+ * layout (file count, structure count, atoms per structure, etc.)
+ *
+ * Example usage:
+ * @code
+ * StructureDB db(...); // fully constructed
+ * auto pinfo = db.exportPass1Info();
+ * // pinfo can be re-used to build a new DB with the same layout
+ * @endcode
+ */
+ inline Pass1Info exportPass1Info() const { return storedPass1_; }
};
} // namespace mlip
diff --git a/include/tadah/mlip/structure_inl.h b/include/tadah/mlip/structure_inl.h
index 55c8a5fadb24bcdabb7d53c2e79e7f19ed095889..2c5e52a6b979782162dd757d63cdfbf9f0889b1e 100644
--- a/include/tadah/mlip/structure_inl.h
+++ b/include/tadah/mlip/structure_inl.h
@@ -61,21 +61,21 @@ inline const double& Structure::fy(size_t i) const { return db_->FY_[offset_ + i
inline const double& Structure::fz(size_t i) const { return db_->FZ_[offset_ + i]; }
inline const int& Structure::atomicNumber(size_t i) const { return db_->atomicNumber_[offset_ + i]; }
-inline Vec3dSoAView Structure::positionView(size_t i)
+inline tadah::core::Vec3dSoAView Structure::positionView(size_t i)
{
- return Vec3dSoAView( x(i), y(i), z(i) );
+ return tadah::core::Vec3dSoAView( x(i), y(i), z(i) );
}
-inline Vec3dSoAView Structure::forceView(size_t i)
+inline tadah::core::Vec3dSoAView Structure::forceView(size_t i)
{
- return Vec3dSoAView( fx(i), fy(i), fz(i) );
+ return tadah::core::Vec3dSoAView( fx(i), fy(i), fz(i) );
}
-inline Vec3dSoAConstView Structure::positionView(size_t i) const
+inline tadah::core::Vec3dSoAConstView Structure::positionView(size_t i) const
{
- return Vec3dSoAConstView( &x(i), &y(i), &z(i) );
+ return tadah::core::Vec3dSoAConstView( &x(i), &y(i), &z(i) );
}
-inline Vec3dSoAConstView Structure::forceView(size_t i) const
+inline tadah::core::Vec3dSoAConstView Structure::forceView(size_t i) const
{
- return Vec3dSoAConstView( &fx(i), &fy(i), &fz(i) );
+ return tadah::core::Vec3dSoAConstView( &fx(i), &fy(i), &fz(i) );
}
inline int Structure::read(std::ifstream &ifs)
@@ -171,7 +171,7 @@ inline double Structure::get_density() const {
V*=1e-24; // convert to cm^3
double amu = 1.66053906660e-24; // grams
double mass = 0;
- for (const auto& a:*this) mass += PeriodicTable::get_mass(a.atomicNumber());
+ for (const auto& a:*this) mass += tadah::core::PeriodicTable::get_mass(a.atomicNumber());
return amu*mass/V;
}
diff --git a/include/tadah/mlip/trainer.h b/include/tadah/mlip/trainer.h
index 8ebc323b3484b1ef55684cc0f5c9d8696a126a19..0396b791769ceb54be3c2fd98b06dfcfa23589db 100644
--- a/include/tadah/mlip/trainer.h
+++ b/include/tadah/mlip/trainer.h
@@ -20,7 +20,7 @@ namespace tadah {
namespace mlip {
class Trainer {
public:
- Config config;
+ tadah::core::Config config;
DC_Selector DCS;
DescriptorsCalc<> dc;
// NNFinder nnf;
@@ -35,15 +35,15 @@ class Trainer {
if(fb)
delete fb;
}
- Trainer (Config &c, tadah::mlip::memory::IMLIPWorkspaceManager& workspaceManager):
+ Trainer (tadah::core::Config &c, tadah::mlip::memory::IMLIPWorkspaceManager& workspaceManager):
config(c),
DCS(config),
dc(config,*DCS.d2b,*DCS.d3b,*DCS.dmb,
*DCS.c2b,*DCS.c3b,*DCS.cmb),
// nnf(config),
- fb(CONFIG::factory<DM_Function_Base,Config&>(
+ fb(tadah::core::factory<DM_Function_Base,tadah::core::Config&>(
config.get<std::string>("MODEL",1),config)),
- model(CONFIG::factory<M_Tadah_Base,DM_Function_Base&,Config&,tadah::mlip::memory::IMLIPWorkspaceManager&>
+ model(tadah::core::factory<M_Tadah_Base,DM_Function_Base&,tadah::core::Config&,tadah::mlip::memory::IMLIPWorkspaceManager&>
(config.get<std::string>("MODEL",0),*fb,config,workspaceManager)),
dm(*fb, config, workspaceManager)
{
@@ -56,7 +56,7 @@ class Trainer {
model->train(stdb,dc);
}
- Config get_param_file() {
+ tadah::core::Config get_param_file() {
return model->get_param_file();
}
@@ -118,7 +118,7 @@ class MPI_Trainer: public Trainer {
int rank;
int ncpu;
- MPI_Trainer(Config &c, int &rank, int &ncpu):
+ MPI_Trainer(tadah::core::Config &c, int &rank, int &ncpu):
Trainer(c),
rank(rank),
ncpu(ncpu)
@@ -291,7 +291,7 @@ class MPI_Trainer: public Trainer {
dm.T.ptr(), &ione, &ione, descB, &context1);
if (rank==0) {
- t_type w(dm.T.ptr(), PHI_cols);
+ tadah::core::t_type w(dm.T.ptr(), PHI_cols);
model->set_weights(w);
// model->trained=true; // still can't train with this model
}
@@ -309,7 +309,7 @@ class TrainerHost: public MPI_Trainer {
std::vector<std::tuple<std::string,int,int>> wpckgs;
public:
- TrainerHost(Config &c, int &rank, int &ncpu):
+ TrainerHost(tadah::core::Config &c, int &rank, int &ncpu):
MPI_Trainer(c, rank, ncpu)
{}
@@ -468,7 +468,7 @@ class TrainerHost: public MPI_Trainer {
class TrainerWorker: public MPI_Trainer {
public:
- TrainerWorker(Config &c, int &rank, int &ncpu):
+ TrainerWorker(tadah::core::Config &c, int &rank, int &ncpu):
MPI_Trainer(c, rank, ncpu)
{}
diff --git a/src/analytics.cpp b/src/analytics.cpp
index d37ac2458feec8b7931968dc28dcadffd8064ed7..48391ea2cd4bc45b3794b237c660580201cfc75d 100644
--- a/src/analytics.cpp
+++ b/src/analytics.cpp
@@ -20,11 +20,11 @@ Analytics::Analytics(const StructureDB &st, const StructureDB &stp)
* @brief Compute per-file MAE of energy (in eV/atom).
* Uses the new approach that sums over all structures in each file f.
*/
-t_type Analytics::calc_e_mae() const
+tadah::core::t_type Analytics::calc_e_mae() const
{
- // We return a t_type with length st.nFiles(),
+ // We return a tadah::core::t_type with length st.nFiles(),
// each element is the MAE for one file.
- t_type emae_vec(st.nFiles());
+ tadah::core::t_type emae_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
@@ -52,9 +52,9 @@ t_type Analytics::calc_e_mae() const
/**
* @brief Compute per-file MAE of forces (in eV/Ã…).
*/
-t_type Analytics::calc_f_mae() const
+tadah::core::t_type Analytics::calc_f_mae() const
{
- t_type fmae_vec(st.nFiles());
+ tadah::core::t_type fmae_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
double fmae = 0.0;
@@ -90,9 +90,9 @@ t_type Analytics::calc_f_mae() const
* @brief Compute per-file MAE of stress (in eV/ų or whichever units).
* We only collect the symmetric components (x,y) for x <= y.
*/
-t_type Analytics::calc_s_mae() const
+tadah::core::t_type Analytics::calc_s_mae() const
{
- t_type smae_vec(st.nFiles());
+ tadah::core::t_type smae_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
double smae = 0.0;
@@ -122,9 +122,9 @@ t_type Analytics::calc_s_mae() const
* @brief Compute per-file relative RMSE of energy:
* sqrt( 1/N * ∑ [ (E1 - E2)/E1 ]² ).
*/
-t_type Analytics::calc_e_rrmse() const
+tadah::core::t_type Analytics::calc_e_rrmse() const
{
- t_type errmse_vec(st.nFiles());
+ tadah::core::t_type errmse_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
double sumSq = 0.0; // sum of squares
@@ -151,9 +151,9 @@ t_type Analytics::calc_e_rrmse() const
/**
* @brief Compute per-file RMSE of energy (in eV/atom).
*/
-t_type Analytics::calc_e_rmse() const
+tadah::core::t_type Analytics::calc_e_rmse() const
{
- t_type ermse_vec(st.nFiles());
+ tadah::core::t_type ermse_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
double sumSq = 0.0;
@@ -179,9 +179,9 @@ t_type Analytics::calc_e_rmse() const
/**
* @brief Compute per-file RMSE of forces (in eV/Ã…).
*/
-t_type Analytics::calc_f_rmse() const
+tadah::core::t_type Analytics::calc_f_rmse() const
{
- t_type frmse_vec(st.nFiles());
+ tadah::core::t_type frmse_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
double sumSq = 0.0;
@@ -212,9 +212,9 @@ t_type Analytics::calc_f_rmse() const
/**
* @brief Compute per-file RMSE of stress (in eV/Ã…^3).
*/
-t_type Analytics::calc_s_rmse() const
+tadah::core::t_type Analytics::calc_s_rmse() const
{
- t_type srmse_vec(st.nFiles());
+ tadah::core::t_type srmse_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
double sumSq = 0.0;
@@ -244,14 +244,14 @@ t_type Analytics::calc_s_rmse() const
/**
* @brief Compute per-file R² of energy (comparing E/atom in each structure).
*/
-t_type Analytics::calc_e_r_sq() const
+tadah::core::t_type Analytics::calc_e_r_sq() const
{
- t_type e_r_sq_vec(st.nFiles());
+ tadah::core::t_type e_r_sq_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
size_t nStructs = st.nStructuresInFile(f);
// Build vectors of observed and predicted (E/atom)
- t_type obs(nStructs), pred(nStructs);
+ tadah::core::t_type obs(nStructs), pred(nStructs);
for (size_t i = 0; i < nStructs; ++i) {
double Eref = st(f,i).energy;
@@ -270,9 +270,9 @@ t_type Analytics::calc_e_r_sq() const
/**
* @brief Compute per-file R² of forces (flattening all atoms’ fx,fy,fz).
*/
-t_type Analytics::calc_f_r_sq() const
+tadah::core::t_type Analytics::calc_f_r_sq() const
{
- t_type f_r_sq_vec(st.nFiles());
+ tadah::core::t_type f_r_sq_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
// First, compute total # of force components (3 * all atoms in file f)
@@ -283,7 +283,7 @@ t_type Analytics::calc_f_r_sq() const
}
size_t nForceComps = 3 * totalAtoms;
- t_type obs(nForceComps), pred(nForceComps);
+ tadah::core::t_type obs(nForceComps), pred(nForceComps);
size_t idx = 0;
// Fill obs/pred with fx, fy, fz
@@ -312,14 +312,14 @@ t_type Analytics::calc_f_r_sq() const
/**
* @brief Compute per-file R² of stress (collecting the 6 independent components).
*/
-t_type Analytics::calc_s_r_sq() const
+tadah::core::t_type Analytics::calc_s_r_sq() const
{
- t_type s_r_sq_vec(st.nFiles());
+ tadah::core::t_type s_r_sq_vec(st.nFiles());
for (size_t f = 0; f < st.nFiles(); ++f)
{
size_t nStructs = st.nStructuresInFile(f);
// Each structure => 6 unique stress components (sxx, syy, szz, sxy, sxz, syz)
- t_type obs(6 * nStructs), pred(6 * nStructs);
+ tadah::core::t_type obs(6 * nStructs), pred(6 * nStructs);
size_t idx = 0;
for (size_t i = 0; i < nStructs; ++i) {
diff --git a/src/castep_castep_reader.cpp b/src/castep_castep_reader.cpp
index 9f3a787ff61304c1704371db686ffecb636e410c..928aa055f962365acc42441687fe735360e6df89 100644
--- a/src/castep_castep_reader.cpp
+++ b/src/castep_castep_reader.cpp
@@ -187,12 +187,12 @@ void CastepCastepReader::parse_data() {
break;
}
- int Z = Element::toAtomicNumber(type);
+ int Z = tadah::core::Element::toAtomicNumber(type);
s->addAtom(Z,px,py,pz,0,0,0);
// s->atoms[i].position = s->cell * s->atoms[i].position; // convert to abs
auto posView = s->positionView(i);
transformBy(posView, s->cell);
- //Vec3d temp = s->cell * s->atoms[i].position; // convert to abs
+ //tadah::core::Vec3d temp = s->cell * s->atoms[i].position; // convert to abs
//s->x(i) = temp(0); s->y(i) = temp(1); s->z(i) = temp(2);
}
if (debug) std::cout << "natoms: " << s->natoms() << std::endl;
@@ -222,7 +222,7 @@ void CastepCastepReader::parse_data() {
std::cerr << "Warning, file" << filename << " line: " << counter << std::endl;
std::cerr << "Warning: Unexpected end of data when reading atom forces: " << line << std::endl;
} else {
- // s->atoms[i].force = Vec3d(fx,fy,fz);
+ // s->atoms[i].force = tadah::core::Vec3d(fx,fy,fz);
s->fx(i) = fx; s->fy(i) = fy; s->fz(i) = fz;
}
}
diff --git a/src/castep_cell_writer.cpp b/src/castep_cell_writer.cpp
index 30255238716f3b2d363598dd79b05e941dbd999d..f8467e8b6db8c44c95cd544e1d1b6a3cd0bed84a 100644
--- a/src/castep_cell_writer.cpp
+++ b/src/castep_cell_writer.cpp
@@ -42,7 +42,7 @@ void CastepCellWriter::write_data(const StructureDB &stdb, const std::string& fi
file << "%BLOCK POSITIONS_ABS" << std::endl;
for (const auto &atom : st) {
- std::string symbol = Element::toSymbol(atom.atomicNumber());
+ std::string symbol = tadah::core::Element::toSymbol(atom.atomicNumber());
file << std::right << std::fixed << std::setw(w)
<< symbol;
file << std::right << std::fixed << std::setw(w)
diff --git a/src/castep_md_reader.cpp b/src/castep_md_reader.cpp
index 92561b2cde23af394cfda30cf812b86d4a53b71c..940082b64148adbc6513d2235ee1b22c01ee3588 100644
--- a/src/castep_md_reader.cpp
+++ b/src/castep_md_reader.cpp
@@ -170,8 +170,8 @@ void CastepMDReader::parse_data() {
std::istringstream iss(line);
if (!(iss >> element >> temp >> px >> py >> pz))
std::cerr << "Warning: Unexpected end of data when reading atomic positions:\n" << line << std::endl;
- // s->add_atom(Atom(Element(element),px*d_conv,py*d_conv,pz*d_conv,0,0,0));
- int Z = Element::toAtomicNumber(element);
+ // s->addAtom(Atom(tadah::core::Element(element),px*d_conv,py*d_conv,pz*d_conv,0,0,0));
+ int Z = tadah::core::Element::toAtomicNumber(element);
s->addAtom(Z,px*d_conv,py*d_conv,pz*d_conv,0,0,0);
R_flag=true;
}
@@ -196,7 +196,7 @@ void CastepMDReader::parse_data() {
F_flag=true;
if (force_idx==s->natoms()) complete_structure=true;
}
- else if (is_blank_line(line)) {
+ else if (tadah::core::is_blank_line(line)) {
if (!error && complete_structure)
postproc_structure(s, structIdxInFile);
diff --git a/src/dm_bf_base.cpp b/src/dm_bf_base.cpp
index ef752f21b7505998204819afc9cd23d45b095f8f..525c8edb5872951a4ec6bf5c13ff47270fbfb1ba 100644
--- a/src/dm_bf_base.cpp
+++ b/src/dm_bf_base.cpp
@@ -4,9 +4,9 @@
namespace tadah {
namespace mlip {
DM_BF_Base::DM_BF_Base() {}
-DM_BF_Base::DM_BF_Base(const Config &c):
- Function_Base(c),
- BF_Base(c),
+DM_BF_Base::DM_BF_Base(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
+ tadah::models::BF_Base(c),
DM_Function_Base(c) {}
DM_BF_Base::~DM_BF_Base() {}
}
diff --git a/src/dm_bf_linear.cpp b/src/dm_bf_linear.cpp
index 43de0a4da9cb1ca22030eef0656ff41be8b03e98..63b3106ecbcd8b3e137bfe9851ad037ad038c355 100644
--- a/src/dm_bf_linear.cpp
+++ b/src/dm_bf_linear.cpp
@@ -3,25 +3,25 @@
namespace tadah {
namespace mlip {
DM_BF_Linear::DM_BF_Linear() {}
-DM_BF_Linear::DM_BF_Linear(const Config &c)
- : Function_Base(c), DM_BF_Base(c), BF_Linear(c) {}
-size_t DM_BF_Linear::get_phi_cols(const Config &config) {
+DM_BF_Linear::DM_BF_Linear(const tadah::core::Config &c)
+ : tadah::models::Function_Base(c), DM_BF_Base(c), tadah::models::BF_Linear(c) {}
+size_t DM_BF_Linear::get_phi_cols(const tadah::core::Config &config) {
size_t cols = config.get<size_t>("DSIZE");
return cols;
}
-void DM_BF_Linear::calc_phi_energy_row(phi_type &Phi, size_t &row,
+void DM_BF_Linear::calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &,
const StructureNeighborView &st_nb,
const StDescriptors &st_d) {
for (size_t a = 0; a < st_d.naed(); ++a) {
- const aed_type &aed = st_d.get_aed(a);
+ const tadah::core::aed_type &aed = st_d.get_aed(a);
for (size_t j = 0; j < aed.size(); ++j) {
Phi(row, j) += aed[j] * fac;
}
}
row++;
}
-void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
+void DM_BF_Linear::calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st,
const StructureNeighborView &st_nb,
const StDescriptors &st_d) {
@@ -29,13 +29,13 @@ void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
for (size_t a=0; a<st.natoms(); ++a) {
const std::size_t* NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- const aed_type& aedi = st_d.get_aed(a);
+ const tadah::core::aed_type& aedi = st_d.get_aed(a);
for (size_t jj=0; jj<NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j= st.globalToLocal(gj); // j is local to st
size_t aa = st_nb.getAAIndex(a, j, jj);
- const fd_type &fij = st_d.fd[a][jj];
- const fd_type &fji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fij = st_d.fd[a][jj];
+ const tadah::core::fd_type &fji = st_d.fd[j][aa];
for (size_t k = 0; k < 3; ++k) {
for (size_t d = 0; d < fij.rows(); ++d) {
Phi(row + k, d) -= fac * (fij(d, k) - fji(d, k));
@@ -45,19 +45,19 @@ void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
row += 3;
}
}
-void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
+void DM_BF_Linear::calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac[6],
const Structure &st,
const StructureNeighborView &st_nb,
const StDescriptors &st_d) {
double V_inv = 1 / st.get_volume();
- Vec3d rj;
+ tadah::core::Vec3d rj;
for (size_t a = 0; a < st.natoms(); ++a) {
const std::size_t *NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- Vec3dSoAConstView ra = st.positionView(a);
- const aed_type &aedi = st_d.get_aed(a);
+ tadah::core::Vec3dSoAConstView ra = st.positionView(a);
+ const tadah::core::aed_type &aedi = st_d.get_aed(a);
for (size_t jj = 0; jj < NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j = st.globalToLocal(gj); // j is local to st
@@ -71,9 +71,9 @@ void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
rj(1) = st.y(j) + sx*st.cell(1,0) + sy*st.cell(1,1) + sz*st.cell(1,2);
rj(2) = st.z(j) + sx*st.cell(2,0) + sy*st.cell(2,1) + sz*st.cell(2,2);
- const fd_type &fdji = st_d.fd[j][aa];
- const fd_type &fdij = st_d.fd[a][jj];
- const aed_type &aedj = st_d.get_aed(j);
+ const tadah::core::fd_type &fdji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fdij = st_d.fd[a][jj];
+ const tadah::core::aed_type &aedj = st_d.get_aed(j);
size_t mn = 0;
for (size_t x = 0; x < 3; ++x) {
diff --git a/src/dm_bf_polynomial2.cpp b/src/dm_bf_polynomial2.cpp
index cf7dbd5cb955359a58f56023f6944ee1b5720088..9bce7b515abea7ce8c3759e9bdcfe7c8695a2eed 100644
--- a/src/dm_bf_polynomial2.cpp
+++ b/src/dm_bf_polynomial2.cpp
@@ -5,24 +5,24 @@
namespace tadah {
namespace mlip {
DM_BF_Polynomial2::DM_BF_Polynomial2() {}
-DM_BF_Polynomial2::DM_BF_Polynomial2(const Config &c):
- Function_Base(c),
+DM_BF_Polynomial2::DM_BF_Polynomial2(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_BF_Base(c),
- BF_Polynomial2(c)
+ tadah::models::BF_Polynomial2(c)
{}
-size_t DM_BF_Polynomial2::get_phi_cols(const Config &config)
+size_t DM_BF_Polynomial2::get_phi_cols(const tadah::core::Config &config)
{
size_t cols = config.get<size_t>("DSIZE");
return (cols*cols+cols)/2;
}
-void DM_BF_Polynomial2::calc_phi_energy_row(phi_type &Phi,
+void DM_BF_Polynomial2::calc_phi_energy_row(tadah::core::phi_type &Phi,
size_t &row,
const double fac,
const Structure &, const StructureNeighborView &st_nb,
const StDescriptors &st_d)
{
for (size_t a=0; a<st_d.naed();++a) {
- const aed_type& aed = st_d.get_aed(a);
+ const tadah::core::aed_type& aed = st_d.get_aed(a);
size_t b=0;
for (size_t i=0; i<st_d.dim(); ++i) {
for (size_t ii=i; ii<st_d.dim(); ++ii) {
@@ -32,7 +32,7 @@ void DM_BF_Polynomial2::calc_phi_energy_row(phi_type &Phi,
}
row++;
}
-void DM_BF_Polynomial2::calc_phi_force_rows(phi_type &Phi,
+void DM_BF_Polynomial2::calc_phi_force_rows(tadah::core::phi_type &Phi,
size_t &row,
const double fac,
const Structure &st, const StructureNeighborView &st_nb,
@@ -42,15 +42,15 @@ void DM_BF_Polynomial2::calc_phi_force_rows(phi_type &Phi,
const std::size_t* NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- Vec3dSoAConstView ra = st.positionView(a);
- const aed_type& aedi = st_d.get_aed(a);
+ tadah::core::Vec3dSoAConstView ra = st.positionView(a);
+ const tadah::core::aed_type& aedi = st_d.get_aed(a);
for (size_t jj=0; jj<NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j= st.globalToLocal(gj); // j is local to st
size_t aa = st_nb.getAAIndex(a, j, jj);
- const fd_type &fdji = st_d.fd[j][aa];
- const fd_type &fdij = st_d.fd[a][jj];
- const aed_type& aedj = st_d.get_aed(j);
+ const tadah::core::fd_type &fdji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fdij = st_d.fd[a][jj];
+ const tadah::core::aed_type& aedj = st_d.get_aed(j);
for (size_t k=0; k<3; ++k) {
size_t b=0;
for (size_t i=0; i<fdij.rows(); ++i) {
@@ -66,20 +66,20 @@ void DM_BF_Polynomial2::calc_phi_force_rows(phi_type &Phi,
row+=3;
}
}
-void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
+void DM_BF_Polynomial2::calc_phi_stress_rows(tadah::core::phi_type &Phi,
size_t &row,
const double fac[6],
const Structure &st, const StructureNeighborView &st_nb,
const StDescriptors &st_d)
{
double V_inv = 1 / st.get_volume();
- Vec3d rj;
+ tadah::core::Vec3d rj;
for (size_t a = 0; a < st.natoms(); ++a) {
const std::size_t *NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- Vec3dSoAConstView ra = st.positionView(a);
- const aed_type &aedi = st_d.get_aed(a);
+ tadah::core::Vec3dSoAConstView ra = st.positionView(a);
+ const tadah::core::aed_type &aedi = st_d.get_aed(a);
for (size_t jj = 0; jj < NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j = st.globalToLocal(gj); // j is local to st
@@ -93,9 +93,9 @@ void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
rj(1) = st.y(j) + sx*st.cell(1,0) + sy*st.cell(1,1) + sz*st.cell(1,2);
rj(2) = st.z(j) + sx*st.cell(2,0) + sy*st.cell(2,1) + sz*st.cell(2,2);
- const fd_type &fdji = st_d.fd[j][aa];
- const fd_type &fdij = st_d.fd[a][jj];
- const aed_type &aedj = st_d.get_aed(j);
+ const tadah::core::fd_type &fdji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fdij = st_d.fd[a][jj];
+ const tadah::core::aed_type &aedj = st_d.get_aed(j);
size_t mn = 0;
for (size_t x = 0; x < 3; ++x) {
for (size_t y = x; y < 3; ++y) {
diff --git a/src/dm_f_all.cpp b/src/dm_f_all.cpp
index 924e3af89b6ed0cbba5695e66431bbc6f2b118cf..43409b220e90c0a0e7b57a4a0a19a900bfe8fa0d 100644
--- a/src/dm_f_all.cpp
+++ b/src/dm_f_all.cpp
@@ -2,25 +2,25 @@
namespace tadah {
namespace mlip {
-template<> CONFIG::Registry<DM_Function_Base>::Map CONFIG::Registry<DM_Function_Base>::registry{};
-template<> CONFIG::Registry<DM_Function_Base,Config&>::Map CONFIG::Registry<DM_Function_Base,Config&>::registry{};
+template<> tadah::core::Registry<DM_Function_Base>::Map tadah::core::Registry<DM_Function_Base>::registry{};
+template<> tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Map tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::registry{};
-CONFIG::Registry<DM_Function_Base>::Register<DM_BF_Linear> DM_BF_Linear_1( "BF_Linear" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_BF_Linear> DM_BF_Linear_2( "BF_Linear" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_BF_Polynomial2> DM_BF_Polynomial2_1( "BF_Polynomial2" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_BF_Polynomial2> DM_BF_Polynomial2_2( "BF_Polynomial2" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Linear> DM_Kern_Linear_1( "Kern_Linear" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Linear> DM_Kern_Linear_2( "Kern_Linear" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_LQ> DM_Kern_LQ_1( "Kern_LQ" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_LQ> DM_Kern_LQ_2( "Kern_LQ" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_1( "Kern_Polynomial" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_2( "Kern_Polynomial" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_1( "Kern_Quadratic" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_2( "Kern_Quadratic" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_RBF> DM_Kern_RBF_1( "Kern_RBF" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_RBF> DM_Kern_RBF_2( "Kern_RBF" );
-CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_1( "Kern_Sigmoid" );
-CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_2( "Kern_Sigmoid" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_BF_Linear> DM_BF_Linear_1( "tadah::models::BF_Linear" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_BF_Linear> DM_BF_Linear_2( "tadah::models::BF_Linear" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_BF_Polynomial2> DM_BF_Polynomial2_1( "tadah::models::BF_Polynomial2" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_BF_Polynomial2> DM_BF_Polynomial2_2( "tadah::models::BF_Polynomial2" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Linear> DM_Kern_Linear_1( "tadah::models::Kern_Linear" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Linear> DM_Kern_Linear_2( "tadah::models::Kern_Linear" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_LQ> DM_Kern_LQ_1( "tadah::models::Kern_LQ" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_LQ> DM_Kern_LQ_2( "tadah::models::Kern_LQ" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_1( "tadah::models::Kern_Polynomial" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_2( "tadah::models::Kern_Polynomial" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_1( "tadah::models::Kern_Quadratic" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_2( "tadah::models::Kern_Quadratic" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_RBF> DM_Kern_RBF_1( "tadah::models::Kern_RBF" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_RBF> DM_Kern_RBF_2( "tadah::models::Kern_RBF" );
+tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_1( "tadah::models::Kern_Sigmoid" );
+tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_2( "tadah::models::Kern_Sigmoid" );
}
}
diff --git a/src/dm_function_base.cpp b/src/dm_function_base.cpp
index 7eb5063d74f666e795f909a28cd4266cab6f4c22..3b408c72f139768a6733f1cef25e7d536808a967 100644
--- a/src/dm_function_base.cpp
+++ b/src/dm_function_base.cpp
@@ -3,7 +3,7 @@
namespace tadah {
namespace mlip {
DM_Function_Base::DM_Function_Base() {}
-DM_Function_Base::DM_Function_Base(const Config &c): Function_Base(c) {}
+DM_Function_Base::DM_Function_Base(const tadah::core::Config &c): tadah::models::Function_Base(c) {}
DM_Function_Base::~DM_Function_Base() {}
}
}
diff --git a/src/dm_kern_base.cpp b/src/dm_kern_base.cpp
index 85d62718e92237a2a2a598824bafb5e703547c42..46baf6034e534e8d02f3c7e0c4590cb8b39023e8 100644
--- a/src/dm_kern_base.cpp
+++ b/src/dm_kern_base.cpp
@@ -7,15 +7,15 @@ namespace tadah {
namespace mlip {
DM_Kern_Base::~DM_Kern_Base() {}
DM_Kern_Base::DM_Kern_Base() {}
-DM_Kern_Base::DM_Kern_Base(const Config &c):
- Function_Base(c),
+DM_Kern_Base::DM_Kern_Base(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
Kern_Base(c),
DM_Function_Base(c) {}
-size_t DM_Kern_Base::get_phi_cols(const Config &)
+size_t DM_Kern_Base::get_phi_cols(const tadah::core::Config &)
{
return basis.cols();
}
-void DM_Kern_Base::calc_phi_energy_row(phi_type &Phi, size_t &row, const double fac,
+void DM_Kern_Base::calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row, const double fac,
const Structure &, const StructureNeighborView &st_nb, const StDescriptors &st_d)
{
for (size_t a=0; a<st_d.naed();++a) {
@@ -26,20 +26,20 @@ void DM_Kern_Base::calc_phi_energy_row(phi_type &Phi, size_t &row, const double
row++;
//Phi.row(row++) *= fac;
}
-void DM_Kern_Base::calc_phi_force_rows(phi_type &Phi, size_t &row, const double fac,
+void DM_Kern_Base::calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row, const double fac,
const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) {
for (size_t a=0; a<st.natoms(); ++a) {
const std::size_t* NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- Vec3dSoAConstView ra = st.positionView(a);
- const aed_type& aedi = st_d.get_aed(a);
+ tadah::core::Vec3dSoAConstView ra = st.positionView(a);
+ const tadah::core::aed_type& aedi = st_d.get_aed(a);
for (size_t jj=0; jj<NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j= st.globalToLocal(gj); // j is local to st
size_t aa = st_nb.getAAIndex(a, j, jj);
- const fd_type &fdji = st_d.fd[j][aa];
- const fd_type &fdij = st_d.fd[a][jj];
- const aed_type& aedj = st_d.get_aed(j);
+ const tadah::core::fd_type &fdji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fdij = st_d.fd[a][jj];
+ const tadah::core::aed_type& aedj = st_d.get_aed(j);
for (size_t b=0; b<basis.cols(); ++b) {
for (size_t k=0; k<3; ++k) {
Phi(row+k,b) -= fac*((*this).prime(basis.col(b), aedi,fdij(k)) -
@@ -50,17 +50,17 @@ void DM_Kern_Base::calc_phi_force_rows(phi_type &Phi, size_t &row, const double
row+=3;
}
}
-void DM_Kern_Base::calc_phi_stress_rows(phi_type &Phi, size_t &row, const double fac[6],
+void DM_Kern_Base::calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row, const double fac[6],
const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d)
{
double V_inv = 1/st.get_volume();
- Vec3d rj;
+ tadah::core::Vec3d rj;
for (size_t a = 0; a < st.natoms(); ++a) {
const std::size_t *NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- Vec3dSoAConstView ra = st.positionView(a);
- const aed_type &aedi = st_d.get_aed(a);
+ tadah::core::Vec3dSoAConstView ra = st.positionView(a);
+ const tadah::core::aed_type &aedi = st_d.get_aed(a);
for (size_t jj = 0; jj < NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j = st.globalToLocal(gj); // j is local to st
@@ -74,9 +74,9 @@ void DM_Kern_Base::calc_phi_stress_rows(phi_type &Phi, size_t &row, const doubl
rj(1) = st.y(j) + sx*st.cell(1,0) + sy*st.cell(1,1) + sz*st.cell(1,2);
rj(2) = st.z(j) + sx*st.cell(2,0) + sy*st.cell(2,1) + sz*st.cell(2,2);
- const fd_type &fdji = st_d.fd[j][aa];
- const fd_type &fdij = st_d.fd[a][jj];
- const aed_type& aedj = st_d.get_aed(j);
+ const tadah::core::fd_type &fdji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fdij = st_d.fd[a][jj];
+ const tadah::core::aed_type& aedj = st_d.get_aed(j);
size_t mn=0;
for (size_t x=0; x<3; ++x) {
for (size_t y=x; y<3; ++y) {
diff --git a/src/dm_kern_linear.cpp b/src/dm_kern_linear.cpp
index 2a675e5b9349074e28409d54953bc92541c3a3d9..67a9020ac5137462457e3c5de7c9bbf8c9450cda 100644
--- a/src/dm_kern_linear.cpp
+++ b/src/dm_kern_linear.cpp
@@ -1,34 +1,34 @@
#include "tadah/mlip/design_matrix/functions/kernels/dm_kern_base.h"
#include <tadah/mlip/design_matrix/functions/kernels/dm_kern_linear.h>
-//CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Linear> DM_Kern_Linear_1( "Kern_Linear" );
-//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Linear> DM_Kern_Linear_2( "Kern_Linear" );
+//tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Linear> DM_Kern_Linear_1( "tadah::models::Kern_Linear" );
+//tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Linear> DM_Kern_Linear_2( "tadah::models::Kern_Linear" );
namespace tadah {
namespace mlip {
DM_Kern_Linear::DM_Kern_Linear() {}
-DM_Kern_Linear::DM_Kern_Linear (const Config &c):
- Function_Base(c),
+DM_Kern_Linear::DM_Kern_Linear (const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_Kern_Base(c),
- Kern_Linear(c)
+ tadah::models::Kern_Linear(c)
{}
-size_t DM_Kern_Linear::get_phi_cols(const Config &config)
+size_t DM_Kern_Linear::get_phi_cols(const tadah::core::Config &config)
{
size_t cols = config.get<size_t>("DSIZE");
return cols;
}
-void DM_Kern_Linear::calc_phi_energy_row(phi_type &Phi, size_t &row,
+void DM_Kern_Linear::calc_phi_energy_row(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &, const StructureNeighborView &st_nb, const StDescriptors &st_d)
{
for (size_t a=0; a<st_d.naed();++a) {
- const aed_type &aed = st_d.get_aed(a); // TODO
+ const tadah::core::aed_type &aed = st_d.get_aed(a); // TODO
for (size_t j=0; j<aed.size(); ++j) {
Phi(row,j)+=aed[j]*fac;
}
}
row++;
}
-void DM_Kern_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
+void DM_Kern_Linear::calc_phi_force_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d)
{
for (size_t a=0; a<st.natoms(); ++a) {
@@ -39,7 +39,7 @@ void DM_Kern_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
const size_t j = st.globalToLocal(gj); // j is local to st
const size_t aa = st_nb.getAAIndex(a, j, jj);
for (size_t k=0; k<3; ++k) {
- aed_type temp = (st_d.fd[a][jj](k)-
+ tadah::core::aed_type temp = (st_d.fd[a][jj](k)-
st_d.fd[j][aa](k))*fac;
for (size_t d=0; d<temp.size(); ++d) {
Phi(row+k,d) -= temp[d];
@@ -50,21 +50,21 @@ void DM_Kern_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
}
}
-void DM_Kern_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
+void DM_Kern_Linear::calc_phi_stress_rows(tadah::core::phi_type &Phi, size_t &row,
const double fac[6], const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d)
{
double V_inv = 1/st.get_volume();
- Vec3d rj;
+ tadah::core::Vec3d rj;
for (size_t a=0; a<st.natoms(); ++a) {
const std::size_t* NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- Vec3dSoAConstView ra = st.positionView(a);
+ tadah::core::Vec3dSoAConstView ra = st.positionView(a);
for (size_t jj=0; jj<NNbrs; ++jj) {
size_t gj = NbPtr[jj];
size_t j= st.globalToLocal(gj); // j is local to st
size_t aa = st_nb.getAAIndex(a, j, jj);
- const fd_type &fdij = st_d.fd[a][jj];
- const fd_type &fdji = st_d.fd[j][aa];
+ const tadah::core::fd_type &fdij = st_d.fd[a][jj];
+ const tadah::core::fd_type &fdji = st_d.fd[j][aa];
auto sx = *st_nb.getShiftXPtr(j);
auto sy = *st_nb.getShiftYPtr(j);
@@ -77,7 +77,7 @@ void DM_Kern_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
size_t mn=0;
for (size_t x=0; x<3; ++x) {
for (size_t y=x; y<3; ++y) {
- aed_type temp = V_inv*(fdij(y)-fdji(y))*0.5*fac[mn]*(ra(x)-rj(x));
+ tadah::core::aed_type temp = V_inv*(fdij(y)-fdji(y))*0.5*fac[mn]*(ra(x)-rj(x));
for (size_t d=0; d<temp.size(); ++d) {
Phi(row+mn,d) += temp[d];
}
diff --git a/src/dm_kern_lq.cpp b/src/dm_kern_lq.cpp
index 5f97ee18933a443e2250780a4ecaa1247684c125..83d42f2f751806508e453e00a54275fe31aa623a 100644
--- a/src/dm_kern_lq.cpp
+++ b/src/dm_kern_lq.cpp
@@ -2,17 +2,17 @@
#include "tadah/models/functions/function_base.h"
#include <tadah/mlip/design_matrix/functions/kernels/dm_kern_lq.h>
-//CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_LQ> DM_Kern_LQ_1( "Kern_LQ" );
-//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_LQ> DM_Kern_LQ_2( "Kern_LQ" );
+//tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_LQ> DM_Kern_LQ_1( "tadah::models::Kern_LQ" );
+//tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_LQ> DM_Kern_LQ_2( "tadah::models::Kern_LQ" );
namespace tadah {
namespace mlip {
DM_Kern_LQ::DM_Kern_LQ()
{}
-DM_Kern_LQ::DM_Kern_LQ(const Config &c):
- Function_Base(c),
+DM_Kern_LQ::DM_Kern_LQ(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_Kern_Base(c),
- Kern_LQ(c)
+ tadah::models::Kern_LQ(c)
{}
}
}
diff --git a/src/dm_kern_polynomial.cpp b/src/dm_kern_polynomial.cpp
index 817597d2d0acf990d23b09dc556dc2973ec04408..f663da0b99f3caf3aaa3227485d128851e3b6594 100644
--- a/src/dm_kern_polynomial.cpp
+++ b/src/dm_kern_polynomial.cpp
@@ -1,16 +1,16 @@
#include <tadah/mlip/design_matrix/functions/kernels/dm_kern_polynomial.h>
-//CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_1( "Kern_Polynomial" );
-//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_2( "Kern_Polynomial" );
+//tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_1( "tadah::models::Kern_Polynomial" );
+//tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Polynomial> DM_Kern_Polynomial_2( "tadah::models::Kern_Polynomial" );
namespace tadah {
namespace mlip {
DM_Kern_Polynomial::DM_Kern_Polynomial()
{}
-DM_Kern_Polynomial::DM_Kern_Polynomial(const Config &c):
- Function_Base(c),
+DM_Kern_Polynomial::DM_Kern_Polynomial(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_Kern_Base(c),
- Kern_Polynomial(c)
+ tadah::models::Kern_Polynomial(c)
{}
}
}
diff --git a/src/dm_kern_quadratic.cpp b/src/dm_kern_quadratic.cpp
index 6ec1aa82d6a53786f2eb06d82461b7a344618809..303ad767679f78690e277c865f4eecbf2076cff1 100644
--- a/src/dm_kern_quadratic.cpp
+++ b/src/dm_kern_quadratic.cpp
@@ -1,17 +1,17 @@
#include "tadah/models/functions/function_base.h"
#include <tadah/mlip/design_matrix/functions/kernels/dm_kern_quadratic.h>
-//CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_1( "Kern_Quadratic" );
-//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_2( "Kern_Quadratic" );
+//tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_1( "tadah::models::Kern_Quadratic" );
+//tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Quadratic> DM_Kern_Quadratic_2( "tadah::models::Kern_Quadratic" );
namespace tadah {
namespace mlip {
DM_Kern_Quadratic::DM_Kern_Quadratic()
{}
-DM_Kern_Quadratic::DM_Kern_Quadratic(const Config &c):
- Function_Base(c),
+DM_Kern_Quadratic::DM_Kern_Quadratic(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_Kern_Base(c),
- Kern_Quadratic(c)
+ tadah::models::Kern_Quadratic(c)
{}
}
}
diff --git a/src/dm_kern_rbf.cpp b/src/dm_kern_rbf.cpp
index 2aa068a7e57e2d94d1384f1fdf3f8a9194cb3ecf..9d2b22505c2f5e8ed4c36c02bf7e169e6e5d38a5 100644
--- a/src/dm_kern_rbf.cpp
+++ b/src/dm_kern_rbf.cpp
@@ -1,16 +1,16 @@
#include <tadah/mlip/design_matrix/functions/kernels/dm_kern_rbf.h>
-//CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_RBF> DM_Kern_RBF_1( "Kern_RBF" );
-//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_RBF> DM_Kern_RBF_2( "Kern_RBF" );
+//tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_RBF> DM_Kern_RBF_1( "tadah::models::Kern_RBF" );
+//tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_RBF> DM_Kern_RBF_2( "tadah::models::Kern_RBF" );
namespace tadah {
namespace mlip {
DM_Kern_RBF::DM_Kern_RBF()
{}
-DM_Kern_RBF::DM_Kern_RBF(const Config &c):
- Function_Base(c),
+DM_Kern_RBF::DM_Kern_RBF(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_Kern_Base(c),
- Kern_RBF(c)
+ tadah::models::Kern_RBF(c)
{}
}
}
diff --git a/src/dm_kern_sigmoid.cpp b/src/dm_kern_sigmoid.cpp
index a75f03109674ffe4f7fc4c13fc9443023f37ff2c..c69250fb8f928d4ea438bb2f437ac49ef63990de 100644
--- a/src/dm_kern_sigmoid.cpp
+++ b/src/dm_kern_sigmoid.cpp
@@ -1,16 +1,16 @@
#include <tadah/mlip/design_matrix/functions/kernels/dm_kern_sigmoid.h>
-//CONFIG::Registry<DM_Function_Base>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_1( "Kern_Sigmoid" );
-//CONFIG::Registry<DM_Function_Base,Config&>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_2( "Kern_Sigmoid" );
+//tadah::core::Registry<DM_Function_Base>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_1( "tadah::models::Kern_Sigmoid" );
+//tadah::core::Registry<DM_Function_Base,tadah::core::Config&>::Register<DM_Kern_Sigmoid> DM_Kern_Sigmoid_2( "tadah::models::Kern_Sigmoid" );
namespace tadah {
namespace mlip {
DM_Kern_Sigmoid::DM_Kern_Sigmoid()
{}
-DM_Kern_Sigmoid::DM_Kern_Sigmoid(const Config &c):
- Function_Base(c),
+DM_Kern_Sigmoid::DM_Kern_Sigmoid(const tadah::core::Config &c):
+ tadah::models::Function_Base(c),
DM_Kern_Base(c),
- Kern_Sigmoid(c)
+ tadah::models::Kern_Sigmoid(c)
{}
}
}
diff --git a/src/lammps_structure_writer.cpp b/src/lammps_structure_writer.cpp
index 86f5468acc3bc30395a9e0a70a73580999ce3f7d..76c87118ed6571bcb8400883752ad997d1269f2d 100644
--- a/src/lammps_structure_writer.cpp
+++ b/src/lammps_structure_writer.cpp
@@ -26,13 +26,13 @@ void LammpsStructureWriter::write_data(const StructureDB &stdb, const std::strin
std::map<std::string, size_t> type;
// Initialize the count for each element
for (const auto& atomicNum : atomicNumbers) {
- std::string symbol = PeriodicTable::find_by_Z(atomicNum).symbol;
+ std::string symbol = tadah::core::PeriodicTable::find_by_Z(atomicNum).symbol;
nelements[symbol] = 0;
}
// then count...
size_t t=0;
for (const auto &atom : st) {
- std::string symbol = PeriodicTable::find_by_Z(atom.atomicNumber()).symbol;
+ std::string symbol = tadah::core::PeriodicTable::find_by_Z(atom.atomicNumber()).symbol;
nelements[symbol]++;
if (type.find(symbol) == type.end()) {
type[symbol] = ++t;
@@ -66,7 +66,7 @@ void LammpsStructureWriter::write_data(const StructureDB &stdb, const std::strin
file << "Masses" << std::endl;
file << std::endl;
for (const auto &t: type) {
- file << t.second << " " << PeriodicTable::get_mass(t.first) << " # " << t.first << std::endl;
+ file << t.second << " " << tadah::core::PeriodicTable::get_mass(t.first) << " # " << t.first << std::endl;
}
file << std::endl;
@@ -76,7 +76,7 @@ void LammpsStructureWriter::write_data(const StructureDB &stdb, const std::strin
size_t idx=1;
for (const auto &atom : st) {
- std::string symbol = PeriodicTable::find_by_Z(atom.atomicNumber()).symbol;
+ std::string symbol = tadah::core::PeriodicTable::find_by_Z(atom.atomicNumber()).symbol;
file << std::right << std::fixed << std::setw(w) << idx;
file << std::right << std::fixed << std::setw(w) << type[symbol];
file << std::right << std::fixed << std::setw(w)
diff --git a/src/m_all.cpp b/src/m_all.cpp
index 5f502149af1710358b2a11532d23e13a31f5a095..ca056d034090190c1a21be29c07bce728d8a0390 100644
--- a/src/m_all.cpp
+++ b/src/m_all.cpp
@@ -5,24 +5,24 @@
namespace tadah {
namespace mlip {
template<>
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&>::Map
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&>::registry{};
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&>::Map
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&>::registry{};
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&>::Register<M_KRR<>> M_KRR_1("M_KRR");
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&>::Register<M_BLR<>> M_BLR_1("M_BLR");
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&>::Register<M_KRR<>> M_KRR_1("M_KRR");
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&>::Register<M_BLR<>> M_BLR_1("M_BLR");
/*template<>*/
-/*CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::models::memory::IModelsWorkspaceManager&>::Map*/
-/*CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::models::memory::IModelsWorkspaceManager&>::registry{};*/
+/*tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::models::memory::IModelsWorkspaceManager&>::Map*/
+/*tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::models::memory::IModelsWorkspaceManager&>::registry{};*/
/**/
-/*CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::models::memory::IModelsWorkspaceManager&>::Register<M_BLR<>> M_BLR_2("M_BLR");*/
-/*CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::models::memory::IModelsWorkspaceManager&>::Register<M_KRR<>> M_KRR_2("M_KRR");*/
+/*tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::models::memory::IModelsWorkspaceManager&>::Register<M_BLR<>> M_BLR_2("M_BLR");*/
+/*tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::models::memory::IModelsWorkspaceManager&>::Register<M_KRR<>> M_KRR_2("M_KRR");*/
template<>
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::Map
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::registry{};
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::Map
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::registry{};
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::Register<M_BLR<>> M_BLR_3("M_BLR");
-CONFIG::Registry<M_Tadah_Base, DM_Function_Base&, Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::Register<M_KRR<>> M_KRR_3("M_KRR");
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::Register<M_BLR<>> M_BLR_3("M_BLR");
+tadah::core::Registry<M_Tadah_Base, DM_Function_Base&, tadah::core::Config&, tadah::mlip::memory::IMLIPWorkspaceManager&>::Register<M_KRR<>> M_KRR_3("M_KRR");
}
}
diff --git a/src/m_tadah_base.cpp b/src/m_tadah_base.cpp
index 51fa45389ac04ec74a7988eabd95b073a982bff7..d53cfe1154ea93a13b380dfc5e1bd0fe50d5fca4 100644
--- a/src/m_tadah_base.cpp
+++ b/src/m_tadah_base.cpp
@@ -9,10 +9,10 @@ fpredict(const size_t a, force_type &v,
for (size_t jj=0; jj<st.nn_size(a); ++jj) {
size_t j=st.near_neigh_idx[a][jj];
const size_t aa = st.get_nn_iindex(a,j,jj);
- const fd_type &fdji = std.fd[j][aa];
- const fd_type &fdij = std.fd[a][jj];
- const aed_type &aedi = std.get_aed(a);
- const aed_type &aedj = std.get_aed(j);
+ const tadah::core::fd_type &fdji = std.fd[j][aa];
+ const tadah::core::fd_type &fdij = std.fd[a][jj];
+ const tadah::core::aed_type &aedi = std.get_aed(a);
+ const tadah::core::aed_type &aedj = std.get_aed(j);
v += fpredict(fdij,aedi);
v -= fpredict(fdji,aedj);
}
@@ -35,7 +35,7 @@ epredict(const StDescriptors &std)
}
Structure M_Tadah_Base::
-predict(const Config &c, StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
+predict(const tadah::core::Config &c, StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
if(c.get<bool>("NORM") && !std.normalised
&& c("MODEL")[1].find("Linear") == std::string::npos)
@@ -61,7 +61,7 @@ predict(const Config &c, StDescriptors &std, const Structure &st, const Structur
}
StructureDB M_Tadah_Base::
-predict(const Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb, const NeighborListDB &nldb)
+predict(const tadah::core::Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb, const NeighborListDB &nldb)
{
StructureDB stdb_;
stdb_.structures.resize(stdb.size());
@@ -76,7 +76,7 @@ predict(const Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb, c
}
StructureDB M_Tadah_Base::
-predict(Config &c, const StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc)
+predict(tadah::core::Config &c, const StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc)
{
StructureDB stdb_;
stdb_.structures.resize(stdb.size());
@@ -114,15 +114,15 @@ spredict(const size_t a, stress_type &s,
const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
double V_inv = 1/st.get_volume();
- const Vec3d &ri = st.atoms[a].position;
- const aed_type &aedi = std.get_aed(a);
+ const tadah::core::Vec3d &ri = st.atoms[a].position;
+ const tadah::core::aed_type &aedi = std.get_aed(a);
for (size_t jj=0; jj<st.nn_size(a); ++jj) {
size_t j=st.near_neigh_idx[a][jj];
const size_t aa = st.get_nn_iindex(a,j,jj);
- const Vec3d &rj = st.nn_pos(a,jj);
- const fd_type &fdij = std.fd[a][jj];
- const fd_type &fdji = std.fd[j][aa];
- const aed_type &aedj = std.get_aed(j);
+ const tadah::core::Vec3d &rj = st.nn_pos(a,jj);
+ const tadah::core::fd_type &fdij = std.fd[a][jj];
+ const tadah::core::fd_type &fdji = std.fd[j][aa];
+ const tadah::core::aed_type &aedj = std.get_aed(j);
const force_type fij = fpredict(fdij,aedi);
const force_type fji = fpredict(fdji,aedj);
@@ -145,15 +145,15 @@ stress_force_predict(const StDescriptors &std, Structure &st_, const StructureNe
double V_inv = 1/st_.get_volume();
for (size_t a=0; a<st_.natoms(); ++a) {
force_type v;
- const Vec3d &ri = st_.atoms[a].position;
- const aed_type &aedi = std.get_aed(a);
+ const tadah::core::Vec3d &ri = st_.atoms[a].position;
+ const tadah::core::aed_type &aedi = std.get_aed(a);
for (size_t jj=0; jj<st_.nn_size(a); ++jj) {
size_t j=st_.near_neigh_idx[a][jj];
const size_t aa = st_.get_nn_iindex(a,j,jj);
- const Vec3d &rj = st_.nn_pos(a,jj);
- const fd_type &fdij = std.fd[a][jj];
- const fd_type &fdji = std.fd[j][aa];
- const aed_type &aedj = std.get_aed(j);
+ const tadah::core::Vec3d &rj = st_.nn_pos(a,jj);
+ const tadah::core::fd_type &fdij = std.fd[a][jj];
+ const tadah::core::fd_type &fdji = std.fd[j][aa];
+ const tadah::core::aed_type &aedj = std.get_aed(j);
const force_type fij = fpredict(fdij,aedi);
const force_type fji = fpredict(fdji,aedj);
diff --git a/src/nn_finder.cpp b/src/nn_finder.cpp
deleted file mode 100644
index 369fa1e399c2145d0c7972294f992e47a7c54103..0000000000000000000000000000000000000000
--- a/src/nn_finder.cpp
+++ /dev/null
@@ -1,290 +0,0 @@
-#include <tadah/mlip/nn_finder.h>
-#include <limits>
-#include <cmath>
-#include <algorithm>
-#include <stdexcept>
-#include <cstring>
-#include <chrono>
-#include <iostream>
-#include <vector>
-
-namespace tadah {
-namespace mlip {
-// Constructor
-NNFinder::NNFinder(Config &config)
-{
- cutoff = config.get<double>("RCUTMAX");
- cutoff_sq = cutoff * cutoff;
-}
-
-// ---------------------------------------------------------------------------
-bool NNFinder::check_box(Structure &st)
-{
- for (int i = 0; i < 3; i++) {
- double vx = st.cell(i, 0);
- double vy = st.cell(i, 1);
- double vz = st.cell(i, 2);
- double len2 = vx*vx + vy*vy + vz*vz;
- if (len2 < cutoff_sq) {
- return false;
- }
- }
- return true;
-}
-
-// ---------------------------------------------------------------------------
-void NNFinder::num_shifts(Structure &st, int N[3]) {
- Matrix3d cell_inv = st.cell.inverse();
-
- double l1 = cell_inv.col(0).norm();
- double l2 = cell_inv.col(1).norm();
- double l3 = cell_inv.col(2).norm();
-
- double f1 = (l1 > 0) ? 1.0/l1 : 1.0;
- double f2 = (l2 > 0) ? 1.0/l2 : 1.0;
- double f3 = (l3 > 0) ? 1.0/l3 : 1.0;
-
- int b1 = std::max(int(f1/cutoff),1);
- int b2 = std::max(int(f2/cutoff),1);
- int b3 = std::max(int(f3/cutoff),1);
-
- N[0] = (int)std::round(0.5 + cutoff*b1/f1);
- N[1] = (int)std::round(0.5 + cutoff*b2/f2);
- N[2] = (int)std::round(0.5 + cutoff*b3/f3);
-}
-
-// ---------------------------------------------------------------------------
-// Naive approach - store only neighbor indices and shifts, no local Atom copies
-void NNFinder::calc_naive(Structure &st)
-{
- st.near_neigh_shift.resize(st.natoms());
- st.near_neigh_idx.resize(st.natoms());
- for (size_t i = 0; i < st.natoms(); i++) {
- st.near_neigh_shift[i].reserve(100);
- st.near_neigh_idx[i].reserve(100);
- }
-
- // Compute shift bounds
- int N[3];
- num_shifts(st, N);
-
- // Precompute shifts (both real displacement and integer triple if needed)
- std::vector<Vec3d> shifts;
- std::vector<Vec3d> shiftIdx;
- shifts.reserve((2*N[0]+1)*(2*N[1]+1)*(2*N[2]+1));
- shiftIdx.reserve((2*N[0]+1)*(2*N[1]+1)*(2*N[2]+1));
-
- for (int n1 = -N[0]; n1 <= N[0]; n1++) {
- for (int n2 = -N[1]; n2 <= N[1]; n2++) {
- for (int n3 = -N[2]; n3 <= N[2]; n3++) {
- Vec3d disp;
- disp[0] = st.cell(0,0)*n1 + st.cell(1,0)*n2 + st.cell(2,0)*n3;
- disp[1] = st.cell(0,1)*n1 + st.cell(1,1)*n2 + st.cell(2,1)*n3;
- disp[2] = st.cell(0,2)*n1 + st.cell(1,2)*n2 + st.cell(2,2)*n3;
- shifts.push_back(disp);
- shiftIdx.push_back(Vec3d(n1, n2, n3));
- }
- }
- }
-
- // Extract positions in contiguous arrays (optional but faster for distance checks)
- const size_t natoms = st.natoms();
- std::vector<double> xPos(natoms), yPos(natoms), zPos(natoms);
- for (size_t i = 0; i < natoms; i++) {
- xPos[i] = st(i).position[0];
- yPos[i] = st(i).position[1];
- zPos[i] = st(i).position[2];
- }
-
- // Distance checks
- for (size_t s = 0; s < shifts.size(); s++) {
- const Vec3d &disp = shifts[s];
- const Vec3d &dispIdxVal = shiftIdx[s];
- bool selfShift = (dispIdxVal[0] == 0 && dispIdxVal[1] == 0 && dispIdxVal[2] == 0);
- size_t startA2 = (selfShift ? 1ul : 0ul);
-
- for (size_t a1 = 0; a1 < natoms; a1++) {
- for (size_t a2 = a1 + startA2; a2 < natoms; a2++) {
- double dx = xPos[a1] - (xPos[a2] + disp[0]);
- double dy = yPos[a1] - (yPos[a2] + disp[1]);
- double dz = zPos[a1] - (zPos[a2] + disp[2]);
- double rij_sq = dx*dx + dy*dy + dz*dz;
- if (rij_sq < cutoff_sq) {
- // forward
- st.near_neigh_idx[a1].push_back(a2);
- st.near_neigh_shift[a1].push_back(dispIdxVal);
-
- // reverse
- st.near_neigh_idx[a2].push_back(a1);
- st.near_neigh_shift[a2].push_back(-dispIdxVal);
- }
- }
- }
- }
-
- // shrink neighbor arrays
- for (size_t i = 0; i < natoms; i++) {
- st.near_neigh_shift[i].shrink_to_fit();
- st.near_neigh_idx[i].shrink_to_fit();
- }
-}
-
-// ---------------------------------------------------------------------------
-// Binned approach - similarly remove local copies, store just indices & shifts
-void NNFinder::calc_binned(Structure &st)
-{
-
- st.near_neigh_shift.resize(st.natoms());
- st.near_neigh_idx.resize(st.natoms());
- for (size_t i = 0; i < st.natoms(); i++) {
- st.near_neigh_shift[i].reserve(100);
- st.near_neigh_idx[i].reserve(100);
- }
-
- // invert cell
- double invC[9];
- inverse_3x3_direct(st.cell.data(), invC);
-
- auto rowLength = [&](int row){
- double vx = st.cell(row, 0);
- double vy = st.cell(row, 1);
- double vz = st.cell(row, 2);
- return std::sqrt(vx*vx + vy*vy + vz*vz);
- };
-
- double cellLenA = rowLength(0);
- double cellLenB = rowLength(1);
- double cellLenC = rowLength(2);
-
- int nBinsA = std::max(1, (int)std::floor(cellLenA / cutoff));
- int nBinsB = std::max(1, (int)std::floor(cellLenB / cutoff));
- int nBinsC = std::max(1, (int)std::floor(cellLenC / cutoff));
-
- struct BinCell {
- std::vector<size_t> atomIndices;
- };
- std::vector<BinCell> bins(nBinsA * nBinsB * nBinsC);
-
- auto binIndex = [&](int ia, int ib, int ic){
- auto wrap = [&](int k, int n){ return ( (k % n) + n ) % n; };
- ia = wrap(ia, nBinsA);
- ib = wrap(ib, nBinsB);
- ic = wrap(ic, nBinsC);
- return size_t(ia + nBinsA * (ib + nBinsB * ic));
- };
-
- // fill bins
- for (size_t i = 0; i < st.natoms(); i++) {
- const auto &atm = st(i);
- double fx = invC[0]*atm.position[0] + invC[3]*atm.position[1] + invC[6]*atm.position[2];
- double fy = invC[1]*atm.position[0] + invC[4]*atm.position[1] + invC[7]*atm.position[2];
- double fz = invC[2]*atm.position[0] + invC[5]*atm.position[1] + invC[8]*atm.position[2];
-
- fx -= std::floor(fx);
- fy -= std::floor(fy);
- fz -= std::floor(fz);
-
- int ia = (int)std::floor(fx*nBinsA);
- int ib = (int)std::floor(fy*nBinsB);
- int ic = (int)std::floor(fz*nBinsC);
- bins[ binIndex(ia, ib, ic) ].atomIndices.push_back(i);
- }
-
- double cutSQ = cutoff_sq;
-
- // search
- for (int ia = 0; ia < nBinsA; ia++) {
- for (int ib = 0; ib < nBinsB; ib++) {
- for (int ic = 0; ic < nBinsC; ic++) {
- size_t b0 = binIndex(ia, ib, ic);
- auto &vec0 = bins[b0].atomIndices;
-
- for (int ja = ia - 1; ja <= ia + 1; ja++) {
- for (int jb = ib - 1; jb <= ib + 1; jb++) {
- for (int jc = ic - 1; jc <= ic + 1; jc++) {
- size_t b1 = binIndex(ja, jb, jc);
- auto &vec1 = bins[b1].atomIndices;
-
- int dA = ja - ia;
- int dB = jb - ib;
- int dC = jc - ic;
-
- for (size_t idxA : vec0) {
- for (size_t idxB : vec1) {
- // avoid double counting
- if ((b0 == b1) && (idxB <= idxA)) continue;
-
- const Atom &a1 = st(idxA);
- const Atom &a2 = st(idxB);
-
- double fracX = double(dA)/double(nBinsA);
- double fracY = double(dB)/double(nBinsB);
- double fracZ = double(dC)/double(nBinsC);
-
- // Real shift
- Vec3d shiftDisp(
- st.cell(0,0)*fracX + st.cell(0,1)*fracY + st.cell(0,2)*fracZ,
- st.cell(1,0)*fracX + st.cell(1,1)*fracY + st.cell(1,2)*fracZ,
- st.cell(2,0)*fracX + st.cell(2,1)*fracY + st.cell(2,2)*fracZ
- );
-
- double dx = a1.position[0] - (a2.position[0] + shiftDisp[0]);
- double dy = a1.position[1] - (a2.position[1] + shiftDisp[1]);
- double dz = a1.position[2] - (a2.position[2] + shiftDisp[2]);
- double dist2 = dx*dx + dy*dy + dz*dz;
-
- if (dist2 < cutSQ) {
- // forward
- st.near_neigh_idx[idxA].push_back(idxB);
- st.near_neigh_shift[idxA].push_back(Vec3d(dA, dB, dC));
-
- // reverse
- st.near_neigh_idx[idxB].push_back(idxA);
- st.near_neigh_shift[idxB].push_back(Vec3d(-dA, -dB, -dC));
- }
- }
- }
- }
- }
- }
- }
- }
- }
-
- // shrink
- for (size_t i = 0; i < st.natoms(); i++) {
- st.near_neigh_shift[i].shrink_to_fit();
- st.near_neigh_idx[i].shrink_to_fit();
- }
-}
-
-// ---------------------------------------------------------------------------
-// Master calc that chooses naive or binned
-void NNFinder::calc(Structure &st)
-{
- if(!check_box(st)) {
- calc_naive(st);
- } else {
- calc_binned(st);
- }
-}
-
-// ---------------------------------------------------------------------------
-// Parallel loop over structure database
-void NNFinder::calc(StructureDB &stdb)
-{
- auto t0 = std::chrono::steady_clock::now();
-#ifdef _OPENMP
- #pragma omp parallel for
-#endif
- for(size_t i = 0; i < stdb.size(); i++){
- calc(stdb(i));
- }
- auto t1 = std::chrono::steady_clock::now();
- double seconds = std::chrono::duration<double>(t1 - t0).count();
- std::cout << "calc(StructureDB &stdb) for-loop took "
- << seconds << " seconds\n";
-}
-
-}
-}
diff --git a/src/st_descriptors.cpp b/src/st_descriptors.cpp
index 25b14f46a439635dc12993179f7dcec829f37627..3378114e292546da2f99089d9942459231012146 100644
--- a/src/st_descriptors.cpp
+++ b/src/st_descriptors.cpp
@@ -2,7 +2,7 @@
namespace tadah {
namespace mlip {
-StDescriptors::StDescriptors(const Structure &s, const Config &c):
+StDescriptors::StDescriptors(const Structure &s, const tadah::core::Config &c):
// fully initialize aed
aeds(c.get<size_t>("DSIZE"), s.natoms()),
// partially init fd
@@ -11,10 +11,10 @@ StDescriptors::StDescriptors(const Structure &s, const Config &c):
//sd(c.get<bool>("STRESS") ? c.get<size_t>("DSIZE"),6 : 0,6)
{
// fd: we still need to resize nn individually for each atom
- // and init fd_type with size
+ // and init tadah::core::fd_type with size
if (c.get<bool>("FORCE") || c.get<bool>("STRESS")) {
for (size_t i=0; i<fd.size(); ++i) {
- fd[i].resize(s.nn_size(i),fd_type(c.get<size_t>("DSIZE")));
+ fd[i].resize(s.nn_size(i),tadah::core::fd_type(c.get<size_t>("DSIZE")));
//for (auto &v:fd[i]) v.set_zero();
}
}
@@ -33,14 +33,14 @@ StDescriptors::StDescriptors(const Structure &s, const Config &c):
}
StDescriptors::StDescriptors() {}
-aed_type & StDescriptors::get_aed(const size_t i) {
+tadah::core::aed_type & StDescriptors::get_aed(const size_t i) {
return aeds.col(i);
}
-const aed_type &StDescriptors::get_aed(const size_t i) const {
+const tadah::core::aed_type &StDescriptors::get_aed(const size_t i) const {
return aeds.col(i);
}
-rho_type& StDescriptors::get_rho(const size_t i) {
+tadah::core::rho_type& StDescriptors::get_rho(const size_t i) {
return rhos.col(i);
}
size_t StDescriptors::naed() const {
diff --git a/src/st_descriptors_db.cpp b/src/st_descriptors_db.cpp
index d7b5b436e9b6b444798c2c114e5ca603bff95045..57c7eed412e4e39b299a97853da711582fea0241 100644
--- a/src/st_descriptors_db.cpp
+++ b/src/st_descriptors_db.cpp
@@ -2,7 +2,7 @@
namespace tadah {
namespace mlip {
-StDescriptorsDB::StDescriptorsDB(const StructureDB &stdb, Config &config):
+StDescriptorsDB::StDescriptorsDB(const StructureDB &stdb, tadah::core::Config &config):
st_descs(stdb.size())
{
// init all StDescriptors
diff --git a/src/vasp_outcar_reader.cpp b/src/vasp_outcar_reader.cpp
index 101b77a49dde7bb67fb1a3540021e21e4ca989d6..df8cec0ad0f67111ff540ea44323f4d5cf6914bd 100644
--- a/src/vasp_outcar_reader.cpp
+++ b/src/vasp_outcar_reader.cpp
@@ -146,7 +146,7 @@ void VaspOutcarReader::parse_data() {
}
size_t count = 0;
for (size_t j = 0; j < natom_types.size(); ++j) {
- Element element(atom_types[j]);
+ tadah::core::Element element(atom_types[j]);
for (size_t i = 0; i < natom_types[j]; ++i) {
if (!std::getline(stream, line)) {
std::cerr << "Warning: Unexpected end of data when reading atom information at atom " << count+1 << std::endl;
@@ -158,7 +158,7 @@ void VaspOutcarReader::parse_data() {
std::cerr << "Warning: Unexpected end of data when reading atom information at atom " << count+1 << std::endl;
} else {
Atom atom(element, px, py, pz, fx, fy, fz);
- s.add_atom(atom);
+ s.addAtom(atom);
count++;
}
}
diff --git a/src/vasp_poscar_writer.cpp b/src/vasp_poscar_writer.cpp
index 42717afb167350b87815c43e2272d2f6521195d7..ca9258046516631f1742d014e5e9ed2b1df2b972 100644
--- a/src/vasp_poscar_writer.cpp
+++ b/src/vasp_poscar_writer.cpp
@@ -11,7 +11,7 @@ VaspPoscarWriter::VaspPoscarWriter() : DatasetWriter() {}
void VaspPoscarWriter::write_data(const StructureDB& stdb, const std::string& filename, size_t structIndex) {
- if (structIndex >= stdb->size()) {
+ if (structIndex >= stdb.size()) {
throw std::out_of_range("Index structIndex is out of range.");
}
diff --git a/src/vasp_vasprun_reader.cpp b/src/vasp_vasprun_reader.cpp
index b8abda8937e02385994ad81fa5f267b880dfbaf2..9b03faed395a87634343ebc534aec4ce67f4d0b4 100644
--- a/src/vasp_vasprun_reader.cpp
+++ b/src/vasp_vasprun_reader.cpp
@@ -201,7 +201,7 @@ void VaspVasprunReader::extract_basis_vectors_and_positions(rx::xml_node<> *stru
if (index < atom_types.size()) {
std::stringstream ss(v_node->value());
if (ss >> x >> y >> z) { // relative positions
- _s.add_atom(Atom(Element(atom_types[index]), x, y, z, 0, 0, 0));
+ _s.addAtom(Atom(tadah::core::Element(atom_types[index]), x, y, z, 0, 0, 0));
_s.atoms.back().position = _s.cell * _s.atoms.back().position; // convert to abs
} else {
std::cerr << "Error parsing atom positions." << std::endl;