diff --git a/include/tadah/mlip/descriptors_calc.hpp b/include/tadah/mlip/descriptors_calc.hpp
index 74b4b91b5b922d12ae501d2ff0b9de27d9fdf884..7cf86c1a321cdc794a57d2d1ef2067b070595509 100644
--- a/include/tadah/mlip/descriptors_calc.hpp
+++ b/include/tadah/mlip/descriptors_calc.hpp
@@ -443,7 +443,9 @@ template <typename D2, typename D3, typename DM, typename C2, typename C3, typen
StDescriptorsDB DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const StructureDB &stdb) {
StDescriptorsDB st_desc_db(stdb, config);
-#pragma omp parallel for
+#ifdef _OPENMP
+ #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));
diff --git a/include/tadah/mlip/design_matrix/design_matrix.h b/include/tadah/mlip/design_matrix/design_matrix.h
index 3f87b5fec3bb6b1188e33da264b93ab0af37c2d2..0289b7a61a3c6ccd522d3bd38419509924968b84 100644
--- a/include/tadah/mlip/design_matrix/design_matrix.h
+++ b/include/tadah/mlip/design_matrix/design_matrix.h
@@ -11,15 +11,15 @@
class DesignMatrixBase {
- public:
- static int phi_rows_num(Config &c, int nstruct_tot, int natoms_tot) {
- bool force = c.get<bool>("FORCE");
- bool stress = c.get<bool>("STRESS");
- int rows = nstruct_tot;
- if (stress) rows += 6*nstruct_tot;
- if (force) rows += 3*natoms_tot;
- return rows;
- }
+public:
+ static int phi_rows_num(Config &c, int nstruct_tot, int natoms_tot) {
+ bool force = c.get<bool>("FORCE");
+ bool stress = c.get<bool>("STRESS");
+ int rows = nstruct_tot;
+ if (stress) rows += 6*nstruct_tot;
+ if (force) rows += 3*natoms_tot;
+ return rows;
+ }
};
@@ -65,19 +65,19 @@ class DesignMatrixBase {
template <typename F>
class DesignMatrix : public DesignMatrixBase {
- public:
+public:
- F f;
- phi_type Phi;
- t_type T;
- t_type Tlabels; // 0-Energy, 1-Force, 2-Stress // TODO should be array of int not double
- bool scale=true; // Control escale,fscale,sscale
+ F f;
+ phi_type Phi;
+ t_type T;
+ t_type Tlabels; // 0-Energy, 1-Force, 2-Stress // TODO should be array of int not double
+ bool scale=true; // Control escale,fscale,sscale
- double e_std_dev=1;
- double f_std_dev=1;
- double s_std_dev[6] = {1,1,1,1,1,1};
+ double e_std_dev=1;
+ double f_std_dev=1;
+ double s_std_dev[6] = {1,1,1,1,1,1};
- /** \brief This constructor fully initialise this object
+ /** \brief This constructor fully initialise this object
*
* This class is used to build a Design Matrix.
*
@@ -89,15 +89,15 @@ class DesignMatrix : public DesignMatrixBase {
*
* \endcode
*/
- DesignMatrix(F &f, Config &c):
- f(f),
- config(c),
- force(config.get<bool>("FORCE")),
- stress(config.get<bool>("STRESS")),
- verbose(config.get<int>("VERBOSE")),
- eweightglob(config.get<double>("EWEIGHT")),
- fweightglob(config.get<double>("FWEIGHT")),
- sweightglob(config.get<double>("SWEIGHT"))
+ DesignMatrix(F &f, Config &c):
+ f(f),
+ config(c),
+ force(config.get<bool>("FORCE")),
+ stress(config.get<bool>("STRESS")),
+ verbose(config.get<int>("VERBOSE")),
+ eweightglob(config.get<double>("EWEIGHT")),
+ fweightglob(config.get<double>("FWEIGHT")),
+ sweightglob(config.get<double>("SWEIGHT"))
{
if (config.exist("ESTDEV"))
e_std_dev = config.get<double>("ESTDEV");
@@ -108,215 +108,219 @@ class DesignMatrix : public DesignMatrixBase {
}
- /** \brief Build design matrix from already calculated StDescriptorsDB
+ /** \brief Build design matrix from already calculated StDescriptorsDB
*
* Here we simply build matrix from already calculated descriptors.
* The vector of targets **T** is build from StructureDB.
* D2, D3 and DM calculators are not used.
*/
- void build(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
- resize(stdb);
- compute_stdevs(stdb);
- fill_T(stdb);
- std::vector<size_t> rows(stdb.size());
- size_t row=0;
- for (size_t s=0; s<stdb.size(); ++s) {
- rows[s] = row;
- row++;
- if (force) {
- for (size_t a=0; a<stdb(s).natoms(); ++a) {
- row+=3;
- }
- }
- if (stress) {
- row+=6;
+ void build(StDescriptorsDB &st_desc_db, const StructureDB &stdb) {
+ resize(stdb);
+ compute_stdevs(stdb);
+ fill_T(stdb);
+ std::vector<size_t> rows(stdb.size());
+ size_t row=0;
+ for (size_t s=0; s<stdb.size(); ++s) {
+ rows[s] = row;
+ row++;
+ if (force) {
+ for (size_t a=0; a<stdb(s).natoms(); ++a) {
+ row+=3;
}
}
- // loop over all structures and build
-#pragma omp parallel for
- for (size_t s=0; s<stdb.size(); ++s) {
- build(rows[s], stdb(s),st_desc_db(s));
+ if (stress) {
+ row+=6;
}
- };
-
- /** \brief Calculate descriptors and build design matrix. */
- template <typename DC>
- void build(const StructureDB &stdb, Normaliser &norm,
- DC &dc) {
- //DescriptorsCalc<D2,D3,DM,C2,C3,CM> dc(config);
- resize(stdb); // call after dc set DSIZE key
- compute_stdevs(stdb);
- fill_T(stdb);
- // for opm we need to find first rows for each structure
-
-
- std::vector<size_t> rows(stdb.size());
- size_t row=0;
- for (size_t s=0; s<stdb.size(); ++s) {
- rows[s] = row;
- row++;
- if (force) {
- for (size_t a=0; a<stdb(s).natoms(); ++a) {
- row+=3;
- }
- }
- if (stress) {
- row+=6;
- }
+ }
+ // loop over all structures and build
+#ifdef _OPENMP
+ #pragma omp parallel for
+#endif
+ for (size_t s=0; s<stdb.size(); ++s) {
+ build(rows[s], stdb(s),st_desc_db(s));
+ }
+ };
+
+ /** \brief Calculate descriptors and build design matrix. */
+ template <typename DC>
+ void build(const StructureDB &stdb, Normaliser &norm,
+ DC &dc) {
+ //DescriptorsCalc<D2,D3,DM,C2,C3,CM> dc(config);
+ resize(stdb); // call after dc set DSIZE key
+ compute_stdevs(stdb);
+ fill_T(stdb);
+ // for opm we need to find first rows for each structure
+
+
+ std::vector<size_t> rows(stdb.size());
+ size_t row=0;
+ for (size_t s=0; s<stdb.size(); ++s) {
+ rows[s] = row;
+ row++;
+ if (force) {
+ for (size_t a=0; a<stdb(s).natoms(); ++a) {
+ row+=3;
}
-#pragma omp parallel for
- for (size_t s=0; s<stdb.size(); ++s) {
- StDescriptors st_d = dc.calc(stdb(s));
+ }
+ if (stress) {
+ row+=6;
+ }
+ }
+#ifdef _OPENMP
+ #pragma omp parallel for
+#endif
+ for (size_t s=0; s<stdb.size(); ++s) {
+ StDescriptors st_d = dc.calc(stdb(s));
- if(config.get<bool>("NORM"))
- norm.normalise(st_d);
+ if(config.get<bool>("NORM"))
+ norm.normalise(st_d);
- build(rows[s], stdb(s), st_d);
- }
+ build(rows[s], stdb(s), st_d);
+ }
- }
- void build(size_t &row, const Structure &st, const StDescriptors &st_d) {
+ }
+ void build(size_t &row, const Structure &st, const StDescriptors &st_d) {
+
+ double escale = 1;
+ if (scale) escale = st.eweight*eweightglob/st.natoms();
+ f.calc_phi_energy_row(Phi,row,escale,st,st_d);
+ if (force) {
+ double fscale = 1;
+ if (scale) fscale = st.fweight*fweightglob/st.natoms()/3.0;///f_std_dev;
+ f.calc_phi_force_rows(Phi,row,fscale,st,st_d);
+ }
+ if (stress) {
+ double sscale_arr[6] {1,1,1,1,1,1};
+ if (scale)
+ for(size_t xy=0;xy<6;++xy)
+ sscale_arr[xy] = st.sweight*sweightglob/6.0;///s_std_dev[xy];
+ f.calc_phi_stress_rows(Phi,row,sscale_arr,st,st_d);
+ }
+ }
+ void fill_T(const StructureDB &stdb, size_t start=0) {
+ size_t j=start;
+ for (size_t s=0; s<stdb.size(); ++s) {
double escale = 1;
- if (scale) escale = st.eweight*eweightglob/st.natoms();
- f.calc_phi_energy_row(Phi,row,escale,st,st_d);
+ if (scale) escale = stdb(s).eweight*eweightglob/stdb(s).natoms();///e_std_dev;
+ Tlabels(j) = 0;
+ T(j++) = stdb(s).energy*escale;
+
if (force) {
double fscale = 1;
- if (scale) fscale = st.fweight*fweightglob/st.natoms()/3.0;///f_std_dev;
- f.calc_phi_force_rows(Phi,row,fscale,st,st_d);
+ if (scale) fscale = stdb(s).fweight*fweightglob/stdb(s).natoms()/3.0;///e_std_dev/f_std_dev;
+ for (const Atom &a : stdb(s).atoms) {
+ Tlabels(j) = 1;
+ T(j++) = a.force(0)*fscale;
+ Tlabels(j) = 1;
+ T(j++) = a.force(1)*fscale;
+ Tlabels(j) = 1;
+ T(j++) = a.force(2)*fscale;
+ }
}
if (stress) {
double sscale_arr[6] {1,1,1,1,1,1};
if (scale)
for(size_t xy=0;xy<6;++xy)
- sscale_arr[xy] = st.sweight*sweightglob/6.0;///s_std_dev[xy];
- f.calc_phi_stress_rows(Phi,row,sscale_arr,st,st_d);
- }
- }
- void fill_T(const StructureDB &stdb, size_t start=0) {
- size_t j=start;
- for (size_t s=0; s<stdb.size(); ++s) {
-
- double escale = 1;
- if (scale) escale = stdb(s).eweight*eweightglob/stdb(s).natoms();///e_std_dev;
- Tlabels(j) = 0;
- T(j++) = stdb(s).energy*escale;
-
- if (force) {
- double fscale = 1;
- if (scale) fscale = stdb(s).fweight*fweightglob/stdb(s).natoms()/3.0;///e_std_dev/f_std_dev;
- for (const Atom &a : stdb(s).atoms) {
- Tlabels(j) = 1;
- T(j++) = a.force(0)*fscale;
- Tlabels(j) = 1;
- T(j++) = a.force(1)*fscale;
- Tlabels(j) = 1;
- T(j++) = a.force(2)*fscale;
- }
- }
- if (stress) {
- double sscale_arr[6] {1,1,1,1,1,1};
- if (scale)
- for(size_t xy=0;xy<6;++xy)
- sscale_arr[xy] = stdb(s).sweight*sweightglob/6.0;///e_std_dev/s_std_dev[xy];
- size_t xy=0;
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- Tlabels(j) = 2;
- T(j++)=stdb(s).stress(x,y)*sscale_arr[xy++];
- }
+ sscale_arr[xy] = stdb(s).sweight*sweightglob/6.0;///e_std_dev/s_std_dev[xy];
+ size_t xy=0;
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ Tlabels(j) = 2;
+ T(j++)=stdb(s).stress(x,y)*sscale_arr[xy++];
}
}
}
}
+ }
- private:
- Config &config;
- size_t rows=0;
- size_t cols=0;
- bool force;
- bool stress;
- int verbose;
- double eweightglob;
- double fweightglob;
- double sweightglob;
-
- // resize Phi and T and set all elements to zero
- void resize(const StructureDB &stdb) {
- rows=0;
- cols = f.get_phi_cols(config);
-
- for (size_t s=0; s<stdb.size(); ++s) {
- rows++;
- if (stress) rows+=6;
- if (force) rows+=stdb(s).natoms()*3;
- }
+private:
+ Config &config;
+ size_t rows=0;
+ size_t cols=0;
+ bool force;
+ bool stress;
+ int verbose;
+ double eweightglob;
+ double fweightglob;
+ double sweightglob;
+
+ // resize Phi and T and set all elements to zero
+ void resize(const StructureDB &stdb) {
+ rows=0;
+ cols = f.get_phi_cols(config);
+
+ for (size_t s=0; s<stdb.size(); ++s) {
+ rows++;
+ if (stress) rows+=6;
+ if (force) rows+=stdb(s).natoms()*3;
+ }
- T.resize(rows);
- Tlabels.resize(rows);
- Phi.resize(rows,cols);
+ T.resize(rows);
+ Tlabels.resize(rows);
+ Phi.resize(rows,cols);
- if (verbose) {
- std::cout << "Phi rows: "<< Phi.rows() << std::endl;
- std::cout << "Phi cols: "<< Phi.cols() << std::endl;
- }
+ if (verbose) {
+ std::cout << "Phi rows: "<< Phi.rows() << std::endl;
+ std::cout << "Phi cols: "<< Phi.cols() << std::endl;
+ }
- T.set_zero();
- Phi.set_zero();
+ T.set_zero();
+ Phi.set_zero();
- }
+ }
- void compute_stdevs(const StructureDB &stdb) {
- t_type evec(stdb.size());
- Matrix svec(stdb.size(),6);
- size_t num_forces=0;
- for (size_t s=0; s<stdb.size(); ++s) {
- evec(s) = stdb(s).energy/stdb(s).natoms();
- if (stress) {
- size_t j=0;
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- svec(s,j++)=stdb(s).stress(x,y);
- }
+ void compute_stdevs(const StructureDB &stdb) {
+ t_type evec(stdb.size());
+ Matrix svec(stdb.size(),6);
+ size_t num_forces=0;
+ for (size_t s=0; s<stdb.size(); ++s) {
+ evec(s) = stdb(s).energy/stdb(s).natoms();
+ if (stress) {
+ size_t j=0;
+ for (size_t x=0; x<3; ++x) {
+ for (size_t y=x; y<3; ++y) {
+ svec(s,j++)=stdb(s).stress(x,y);
}
}
- num_forces +=3*stdb(s).natoms();
}
+ num_forces +=3*stdb(s).natoms();
+ }
- //e_std_dev = std::sqrt((evec - evec.mean()).square().sum()/(evec.size()-1));
- e_std_dev = evec.std_dev(evec.mean(), evec.size()-1);
- if (verbose) std::cout << "Energy standard deviation (eV/atom): " << e_std_dev << std::endl;
+ //e_std_dev = std::sqrt((evec - evec.mean()).square().sum()/(evec.size()-1));
+ e_std_dev = evec.std_dev(evec.mean(), evec.size()-1);
+ if (verbose) std::cout << "Energy standard deviation (eV/atom): " << e_std_dev << std::endl;
- if (stress) {
- for (size_t j=0; j<6; ++j) {
- s_std_dev[j] = svec.col(j).std_dev(svec.col(j).mean(), svec.col(j).size()-1);
- if (verbose) std::cout << "Stress standard deviation (eV): " << s_std_dev[j] << std::endl;
- }
+ if (stress) {
+ for (size_t j=0; j<6; ++j) {
+ s_std_dev[j] = svec.col(j).std_dev(svec.col(j).mean(), svec.col(j).size()-1);
+ if (verbose) std::cout << "Stress standard deviation (eV): " << s_std_dev[j] << std::endl;
}
+ }
- if (force) {
- size_t j=0;
- t_type fvec(num_forces);
- for (size_t s=0; s<stdb.size(); ++s) {
- for (const Atom &a : stdb(s).atoms) {
- fvec(j++) = a.force(0);
- fvec(j++) = a.force(1);
- fvec(j++) = a.force(2);
- }
+ if (force) {
+ size_t j=0;
+ t_type fvec(num_forces);
+ for (size_t s=0; s<stdb.size(); ++s) {
+ for (const Atom &a : stdb(s).atoms) {
+ fvec(j++) = a.force(0);
+ fvec(j++) = a.force(1);
+ fvec(j++) = a.force(2);
}
- //fvec /= e_std_dev;
- //svec /= e_std_dev;
- // e_std_dev has units of energy
- // f_std_dev has units of inverse distance
- f_std_dev = fvec.std_dev(fvec.mean(),fvec.size()-1);
- if (verbose) std::cout << "Force standard deviation (1/A): " << f_std_dev << std::endl;
}
-
- config.add("ESTDEV",e_std_dev);
- config.add("FSTDEV",f_std_dev);
- for (size_t j=0; j<6; ++j)
- config.add("SSTDEV",s_std_dev[j]);
+ //fvec /= e_std_dev;
+ //svec /= e_std_dev;
+ // e_std_dev has units of energy
+ // f_std_dev has units of inverse distance
+ f_std_dev = fvec.std_dev(fvec.mean(),fvec.size()-1);
+ if (verbose) std::cout << "Force standard deviation (1/A): " << f_std_dev << std::endl;
}
+
+ config.add("ESTDEV",e_std_dev);
+ config.add("FSTDEV",f_std_dev);
+ for (size_t j=0; j<6; ++j)
+ config.add("SSTDEV",s_std_dev[j]);
+ }
};
#endif
diff --git a/include/tadah/mlip/structure_db.h b/include/tadah/mlip/structure_db.h
index 74697d6edb2dc255ab051141ea135e12337348b6..35c8f197c798e2fd38b436958a0d85d8e6d84d31 100644
--- a/include/tadah/mlip/structure_db.h
+++ b/include/tadah/mlip/structure_db.h
@@ -73,7 +73,7 @@ struct StructureDB {
*
* Return the number of loaded structures.
*/
- int add(const std::string fn, int first, int N);
+ int add(const std::string fn, size_t first, int N);
/** Add a single Structure object to this container */
void add(const Structure &s);
diff --git a/src/m_tadah_base.cpp b/src/m_tadah_base.cpp
index d78f49839c6cea61e2c57595059d130f0bfb280a..66a70b0747838608dceb6505fbefeb4caa328e37 100644
--- a/src/m_tadah_base.cpp
+++ b/src/m_tadah_base.cpp
@@ -63,7 +63,9 @@ predict(const Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb)
{
StructureDB stdb_;
stdb_.structures.resize(stdb.size());
+#ifdef _OPENMP
#pragma omp parallel for
+#endif
for (size_t i=0; i<stdb.size(); ++i) {
stdb_(i) = predict(c,st_desc_db(i), stdb(i));
}
@@ -75,7 +77,9 @@ predict(Config &c, const StructureDB &stdb, DC_Base &dc)
{
StructureDB stdb_;
stdb_.structures.resize(stdb.size());
+#ifdef _OPENMP
#pragma omp parallel for
+#endif
for (size_t i=0; i<stdb.size(); ++i) {
StDescriptors st_d = dc.calc(stdb(i));
stdb_(i) = predict(c,st_d, stdb(i));
diff --git a/src/nn_finder.cpp b/src/nn_finder.cpp
index 9c41cb4d4fbe4713447c104111ecbf7fc5e8bfdf..c2840d9cd04d2265724e183075f772b6d8f2819c 100644
--- a/src/nn_finder.cpp
+++ b/src/nn_finder.cpp
@@ -72,7 +72,9 @@ void NNFinder::calc(Structure &st) {
}
}
void NNFinder::calc(StructureDB &stdb) {
+#ifdef _OPENMP
#pragma omp parallel for
+#endif
for (size_t s=0; s<stdb.size(); ++s) {
calc(stdb(s));
}
diff --git a/src/structure_db.cpp b/src/structure_db.cpp
index 08caaffbb156c3b4967c53675ab1515c3822fbf0..5966f683c35612380528a0d9f1e1a2835aca9e80 100644
--- a/src/structure_db.cpp
+++ b/src/structure_db.cpp
@@ -27,7 +27,7 @@ void StructureDB::add(const std::string fn) {
}
ifs.close();
}
-int StructureDB::add(const std::string fn, int first, int N) {
+int StructureDB::add(const std::string fn, size_t first, int N) {
std::ifstream ifs(fn);
if (!ifs.is_open()) {
throw std::runtime_error("DBFILE does not exist: "+fn);