diff --git a/include/tadah/mlip/neighbor_list_db.h b/include/tadah/mlip/neighbor_list_db.h
index acd4ea43ddce7a012f84f0cd619d9c70ed576081..9c062e119aae85f62315bd3ae08e6abe9a8eaccb 100644
--- a/include/tadah/mlip/neighbor_list_db.h
+++ b/include/tadah/mlip/neighbor_list_db.h
@@ -5,7 +5,6 @@
#include <tadah/mlip/neighbor_list.h>
#include <tadah/mlip/structure_db.h>
-// #include <tadah/mlip/structure_neighbor_view.h>
namespace tadah {
namespace mlip {
diff --git a/include/tadah/mlip/structure.h b/include/tadah/mlip/structure.h
index 5bafff0c7cfba27b8f00edc1f41aab1a1fb7d1a5..3ebe47f34f8730cb83b38c8ce066d044d0f591e7 100644
--- a/include/tadah/mlip/structure.h
+++ b/include/tadah/mlip/structure.h
@@ -12,17 +12,17 @@
#include <tadah/core/periodic_table.h>
#include <tadah/core/core_types.h>
#include <tadah/core/vec3d_soa_view.h>
+#include <tadah/core/vec3d_soa_const_view.h>
+
+// TODO This class must be reworked such that it does not store any data
+// It should be a lightweight handle similar to how Atom works.
-// Forward-declare classes to avoid circular includes
namespace tadah {
namespace mlip {
+
+// Forward-declare classes to avoid circular includes
class Atom; // We'll define it in "atom.h"
class StructureDB; // We'll define it in "structure_db.h"
-} // namespace mlip
-} // namespace tadah
-
-namespace tadah {
-namespace mlip {
/**
* @brief Represents a subrange of SoA arrays for storing atoms.
@@ -87,6 +87,9 @@ public:
Vec3dSoAView positionView(size_t i);
Vec3dSoAView forceView(size_t i);
+ Vec3dSoAConstView positionView(size_t i) const;
+ Vec3dSoAConstView forceView(size_t i) const;
+
/**
* @brief Converts a global index to a local index within this structure.
*
@@ -219,4 +222,3 @@ private:
} // namespace tadah
#endif // TADAH_MLIP_STRUCTURE_H
-
diff --git a/include/tadah/mlip/structure_inl.h b/include/tadah/mlip/structure_inl.h
index b49cfb89dc4b02e3ccb97ee27e4dd73cf226062a..55c8a5fadb24bcdabb7d53c2e79e7f19ed095889 100644
--- a/include/tadah/mlip/structure_inl.h
+++ b/include/tadah/mlip/structure_inl.h
@@ -69,6 +69,14 @@ inline Vec3dSoAView Structure::forceView(size_t i)
{
return Vec3dSoAView( fx(i), fy(i), fz(i) );
}
+inline Vec3dSoAConstView Structure::positionView(size_t i) const
+{
+ return Vec3dSoAConstView( &x(i), &y(i), &z(i) );
+}
+inline Vec3dSoAConstView Structure::forceView(size_t i) const
+{
+ return Vec3dSoAConstView( &fx(i), &fy(i), &fz(i) );
+}
inline int Structure::read(std::ifstream &ifs)
{
diff --git a/include/tadah/mlip/structure_neighbor_view.h b/include/tadah/mlip/structure_neighbor_view.h
index 21b9b1b11d684e849dc72006ff6713276477acbb..fb871478d2b54c0ab262d3475db2210b5cd7d0b7 100644
--- a/include/tadah/mlip/structure_neighbor_view.h
+++ b/include/tadah/mlip/structure_neighbor_view.h
@@ -91,7 +91,7 @@ public:
}
/**
- * @brief Returns the X-shift pointer, if HPC code uses shifts.
+ * @brief Returns the X-shift pointer.
*/
inline const int* getShiftXPtr(std::size_t localAtomIndex) const {
if (localAtomIndex >= str_.natoms()) {
diff --git a/src/dm_bf_linear.cpp b/src/dm_bf_linear.cpp
index 483c263d7c330f7f25c9ccec1e835406f854bf4f..77f84db27b2015bdc5aed340ec33ebab4c2a7f7c 100644
--- a/src/dm_bf_linear.cpp
+++ b/src/dm_bf_linear.cpp
@@ -1,68 +1,75 @@
#include <tadah/mlip/design_matrix/functions/basis_functions/dm_bf_linear.h>
-//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" );
-
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 Config &c)
+ : Function_Base(c), DM_BF_Base(c), BF_Linear(c) {}
+size_t DM_BF_Linear::get_phi_cols(const Config &config) {
size_t cols = config.get<size_t>("DSIZE");
return cols;
}
void DM_BF_Linear::calc_phi_energy_row(phi_type &Phi, size_t &row,
- const double fac, const Structure &, const StructureNeighborView &st_nb, const StDescriptors &st_d)
-{
- for (size_t i=0; i<st_d.naed(); ++i) {
- const aed_type &aed = st_d.get_aed(i);
- for (size_t j=0; j<aed.size(); ++j) {
- Phi(row,j)+=aed[j]*fac;
+ 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);
+ for (size_t j = 0; j < aed.size(); ++j) {
+ Phi(row, j) += aed[j] * fac;
}
}
row++;
}
void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
- const double fac, const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d)
-{
+ const double fac, const Structure &st,
+ const StructureNeighborView &st_nb,
+ const StDescriptors &st_d) {
- for (size_t a=0; a<st.natoms(); ++a) {
- for (size_t jj=0; jj<st_d.fd[a].size(); ++jj) {
- const size_t j=st.near_neigh_idx[a][jj];
- const size_t aa = st.get_nn_iindex(a,j,jj);
+ 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 Vec3dSoAConstView ri = st.positionView(a);
+ const 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 = getAAIndex(a, j, jj)
const fd_type &fij = st_d.fd[a][jj];
const fd_type &fji = st_d.fd[j][aa];
- for (size_t k=0; k<3; ++k) {
- for (size_t d=0; d<fij.rows(); ++d) {
- Phi(row+k,d) -= fac*(fij(d,k)-fji(d,k));
+ for (size_t k = 0; k < 3; ++k) {
+ for (size_t d = 0; d < fij.rows(); ++d) {
+ Phi(row + k, d) -= fac * (fij(d, k) - fji(d, k));
}
}
}
- row+=3;
+ row += 3;
}
}
void DM_BF_Linear::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)
-{
- double V_inv = 1/st.get_volume();
- for (size_t i=0; i<st.natoms(); ++i) {
- const Vec3d &ri = st(i).position;
- for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
- const size_t j=st.near_neigh_idx[i][jj];
- const size_t ii = st.get_nn_iindex(i,j,jj);
- const fd_type &fdij = st_d.fd[i][jj];
- const fd_type &fdji = st_d.fd[j][ii];
- const Vec3d &rj = st.nn_pos(i,jj);
- size_t mn=0;
- for (size_t x=0; x<3; ++x) {
- for (size_t y=x; y<3; ++y) {
- for (size_t d=0; d<fdij.rows(); ++d) {
- Phi(row+mn,d) += V_inv*(fdij(d,y)-fdji(d,y))*0.5*fac[mn]*(ri(x)-rj(x));
+ const double fac[6],
+ const Structure &st,
+ const StructureNeighborView &st_nb,
+ const StDescriptors &st_d) {
+ double V_inv = 1 / st.get_volume();
+ 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 Vec3dSoAConstView ri = st.positionView(a);
+ const 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 = getAAIndex(a, j, jj)
+ const fd_type &fdij = st_d.fd[a][jj];
+ const fd_type &fdji = st_d.fd[j][aa];
+ const Vec3d &rj = st.nn_pos(a, jj);
+ size_t mn = 0;
+ for (size_t x = 0; x < 3; ++x) {
+ for (size_t y = x; y < 3; ++y) {
+ for (size_t d = 0; d < fdij.rows(); ++d) {
+ Phi(row + mn, d) += V_inv * (fdij(d, y) - fdji(d, y)) * 0.5 *
+ fac[mn] * (ri(x) - rj(x));
}
mn++;
}
@@ -71,5 +78,5 @@ void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
}
row += 6;
}
-}
-}
+} // namespace mlip
+} // namespace tadah
diff --git a/src/dm_bf_polynomial2.cpp b/src/dm_bf_polynomial2.cpp
index 24f22e31975963a5cffd4d9304d47b6ffe7570b2..b139f79e51b172244ae9e8d67feff18a8b468c2c 100644
--- a/src/dm_bf_polynomial2.cpp
+++ b/src/dm_bf_polynomial2.cpp
@@ -2,9 +2,6 @@
#include "tadah/models/functions/function_base.h"
#include <tadah/mlip/design_matrix/functions/basis_functions/dm_bf_polynomial2.h>
-//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" );
-
namespace tadah {
namespace mlip {
DM_BF_Polynomial2::DM_BF_Polynomial2() {}
@@ -45,7 +42,7 @@ 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);
- const Vec3dSoAView ri = st.positionView(a);
+ const Vec3dSoAView ra = st.positionView(a);
const aed_type& aedi = st_d.get_aed(a);
for (size_t jj=0; jj<NNbrs; ++jj) {
@@ -79,32 +76,40 @@ void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
const Structure &st, const StructureNeighborView &st_nb,
const StDescriptors &st_d)
{
- double V_inv = 1/st.get_volume();
- for (size_t a=0; a<st.natoms(); ++a) {
+ double V_inv = 1 / st.get_volume();
+ for (size_t a = 0; a < st.natoms(); ++a) {
- const std::size_t* NbPtr = st_nb.getNeighborsPtr(a);
+ const std::size_t *NbPtr = st_nb.getNeighborsPtr(a);
size_t NNbrs = st_nb.numNeighbors(a);
- const Vec3dSoAView ri = st.positionView(a);
- const aed_type& aedi = st_d.get_aed(a);
+ const Vec3dSoAView ra = st.positionView(a);
+ const aed_type &aedi = st_d.get_aed(a);
+ Vec3d rj;
+ 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 = getAAIndex(a, j, jj);
- for (size_t jj=0; jj<NNbrs; ++jj) {
+ auto sx = st.getShiftXPtr(j);
+ auto sy = st.getShiftYPtr(j);
+ auto sz = st.getShiftZPtr(j);
+
+ rj(0) = st(j).x() + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ rj(1) = st(j).y() + sx*st.cell(1,0) + sy*st.cell(1,1) + sz*st.cell(1,2);
+ rj(2) = st(j).z() + sx*st.cell(2,0) + sy*st.cell(2,1) + sz*st.cell(2,2);
- size_t gj = NbPtr[jj];
- size_t j= st.globalToLocal(gj); // j is local to st
- size_t aa = getAAIndex(a, j, jj)
const fd_type &fdji = st_d.fd[j][aa];
const fd_type &fdij = st_d.fd[a][jj];
- const Vec3d &rj = st.nn_pos(a,jj);
- const aed_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) {
- size_t b=0;
- for (size_t i=0; i<fdij.rows(); ++i) {
- for (size_t ii=i; ii<fdij.rows(); ++ii) {
- Phi(row+mn,b) += V_inv*0.5*fac[mn]*(ri(x)-rj(x))
- *(fdij(i,y)*aedi(ii) + fdij(ii,y)*aedi(i)
- - fdji(i,y)*aedj(ii) - fdji(ii,y)*aedj(i));
+ const 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) {
+ size_t b = 0;
+ for (size_t i = 0; i < fdij.rows(); ++i) {
+ for (size_t ii = i; ii < fdij.rows(); ++ii) {
+ Phi(row + mn, b) +=
+ V_inv * 0.5 * fac[mn] * (ra(x) - rj(x)) *
+ (fdij(i, y) * aedi(ii) + fdij(ii, y) * aedi(i) -
+ fdji(i, y) * aedj(ii) - fdji(ii, y) * aedj(i));
b++;
}
}
diff --git a/src/structure.cpp b/src/structure.cpp
deleted file mode 100644
index ee3ca19c640ae715b5e7d92070e5cebadf3e023d..0000000000000000000000000000000000000000
--- a/src/structure.cpp
+++ /dev/null
@@ -1,354 +0,0 @@
-#include <tadah/mlip/structure.h>
-#include <tadah/core/periodic_table.h>
-
-#include <stdexcept>
-#include <cmath>
-
-Structure::Structure() {
- PeriodicTable::initialize();
-}
-
-const Atom &Structure::operator()(const size_t i) const{
- return atoms[i];
-}
-
-void Structure::add_atom(Atom a) {
- atoms.push_back(a);
-}
-void Structure::remove_atom(const size_t i) {
- atoms.erase(atoms.begin()+i);
-}
-
-Vec3d Structure::nn_pos(const size_t i, const size_t n) const
-{
- // (A) Global index of this neighbor
- const size_t neighborIndex = near_neigh_idx[i][n];
-
- // (B) Atom's original "unshifted" position:
- const Vec3d &posNeighbor = atoms[neighborIndex].position;
-
- // (C) Convert the stored shift -> real displacement shiftDisp
- // If near_neigh_shift[i][n] is an integer triple (n1, n2, n3),
- // multiply by the cell. Otherwise, if it's already in real space,
- // you can just do: Vec3d shiftDisp = near_neigh_shift[i][n].
- Vec3d shift = near_neigh_shift[i][n]; // might be (n1, n2, n3)
-
- // If shift is integer triple (n1, n2, n3), multiply with cell:
- Vec3d shiftDisp;
- shiftDisp[0] = shift[0]*cell(0,0) + shift[1]*cell(0,1) + shift[2]*cell(0,2);
- shiftDisp[1] = shift[0]*cell(1,0) + shift[1]*cell(1,1) + shift[2]*cell(1,2);
- shiftDisp[2] = shift[0]*cell(2,0) + shift[1]*cell(2,1) + shift[2]*cell(2,2);
-
- // (D) Final neighbor position = unshifted position + shiftDisp
- Vec3d posShifted = posNeighbor + shiftDisp;
- return posShifted;
-}
-
-size_t Structure::natoms() const {
- return atoms.size();
-}
-
-size_t Structure::nn_size(size_t i) const {
- return near_neigh_idx[i].size();
-}
-
-int Structure::read(std::ifstream &ifs) {
- std::string line;
- // ignore extra empty lines
- std::getline(ifs,line);
- if(line.empty()) return 1;
-
- // OK, process structural data
- // read first line as a label
- label = line;
-
- // the second line could be energy or
- // a scalling factors eweight fweight sweight
- std::getline(ifs,line);
- std::stringstream stream(line);
- std::string temp;
- size_t count = 0;
- while(stream >> temp) { ++count;}
-
- if (count == 3) {
- stream.clear();
- stream.seekg(0, std::ios::beg);
- stream >> eweight >> fweight >> sweight;
- std::getline(ifs,line);
- std::stringstream stream2(line);
- // energy and T
- stream2 >> energy;
- if(stream2) stream2 >> T;
- }
- else if (count == 2) {
- stream.clear();
- stream.seekg(0, std::ios::beg);
- // energy and T
- stream >> energy >> T;
- }
- else {
- energy = std::stod(line);
- }
-
- // 3 lines, 9 elements are for the cell matrix
- for (int i=0; i<3; ++i)
- for (int j=0; j<3; ++j)
- ifs >> cell(i,j);
-
- // 3 lines, 9 elements are for the stress matrix
- // if (use_stress)
- for (int i=0; i<3; ++i)
- for (int j=0; j<3; ++j)
- ifs >> stress(i,j);
-
- // move to next line
- std::getline(ifs,line);
-
- // make sure atoms vector is empty
- atoms.clear();
-
- char symbol[3];
- double px,py,pz,fx,fy,fz;
- while(std::getline(ifs,line)) {
- if(line.empty()) break;
- //if(line == " ") break;
- if(line == "\r") break; // detects windows newline
-
- std::istringstream tmp(line);
- tmp >> symbol >> px >> py >> pz >> fx >> fy >> fz;
- Element e = PeriodicTable::find_by_symbol(symbol);
- atoms.push_back(Atom(e,px,py,pz,fx,fy,fz));
- //Atom &a = atoms.back();
-
- //tmp >> a.label >> a.position(0) >> a.position(1) >> a.position(2);
- // if (use_force)
- //tmp >> a.force(0) >> a.force(1) >> a.force(2);
- //std::cout << a.force(0) << " " << a.force(1) << " " << a.force(2) << std::endl;;
- }
- return 0;
-}
-
-void Structure::read(const std::string fn) {
- std::ifstream ifs(fn);
- if (!ifs.is_open()) {
- throw std::runtime_error("File does not exist: "+fn);
- }
- read(ifs);
- ifs.close();
-}
-
-void Structure::save(std::ofstream &ofs) const {
- ofs << label << std::endl;
- ofs << eweight << " " << fweight << " " << sweight << std::endl;
- ofs << energy << std::endl;
- for (int i=0; i<3; ++i) {
- ofs << std::fixed << std::right << std::setw(w) << std::setprecision(p) << cell(i,0);
- ofs << std::fixed << std::right << std::setw(w) << std::setprecision(p) << cell(i,1);
- ofs << std::fixed << std::right << std::setw(w) << std::setprecision(p) << cell(i,2);
- ofs << std::endl;
- }
-
- for (int i=0; i<3; ++i) {
- ofs << std::fixed << std::right << std::setw(w) << std::setprecision(p) << stress(i,0);
- ofs << std::fixed << std::right << std::setw(w) << std::setprecision(p) << stress(i,1);
- ofs << std::fixed << std::right << std::setw(w) << std::setprecision(p) << stress(i,2);
- ofs << std::endl;
- }
-
- for (const Atom &a:atoms) {
- ofs << a.symbol[0] << a.symbol[1];
- ofs << std::setw(w-2) << std::fixed << std::right
- << std::setprecision(p) << a.position(0);
- ofs << std::setw(w) << std::fixed << std::right
- << std::setprecision(p) << a.position(1);
- ofs << std::setw(w) << std::fixed << std::right
- << std::setprecision(p) << a.position(2);
- ofs << std::setw(w) << std::fixed << std::right
- << std::setprecision(p) << a.force(0);
- ofs << std::setw(w) << std::fixed << std::right
- << std::setprecision(p) << a.force(1);
- ofs << std::setw(w) << std::fixed << std::right
- << std::setprecision(p) << a.force(1);
- ofs << std::endl;
- }
-}
-
-void Structure::save(const std::string fn) const {
- std::ofstream ofs(fn);
- save(ofs);
- ofs.close();
-}
-size_t Structure::get_nn_iindex(const size_t i, const size_t j, const size_t jj) const {
- Vec3d shift_ijj = -near_neigh_shift[i][jj];
- size_t ii = 0;
-
- for (size_t x=0; x<near_neigh_idx[j].size(); ++x) {
- if (near_neigh_idx[j][x] == i)
- if (near_neigh_shift[j][ii] == shift_ijj)
- break;
- ii++;
- }
- return ii;
-
-}
-double Structure::get_volume() const {
- return cell.row(0)*(cell.row(1).cross(cell.row(2)));
-}
-double Structure::get_density() const {
- double V = get_volume();
- V*=1e-24; // convert to cm^3
- double amu = 1.66053906660e-24; // g
- double mass = 0;
- for (const auto& a:atoms) mass += PeriodicTable::get_mass(a.Z);
- return amu*mass/V;
-}
-
-double Structure::get_temperature() const {
- return T;
-}
-
-double Structure::get_virial_pressure() const {
- return stress.trace()/get_volume()/3;
-}
-
-double Structure::get_pressure(double T, double kB) const {
- double vpress = get_virial_pressure();
- return vpress + natoms()*kB*T/get_volume();
-}
-bool Structure::operator==(const Structure& st) const
-{
- double EPSILON = std::numeric_limits<double>::epsilon();
- bool result =
- cell.isApprox(st.cell, EPSILON)
- && stress.isApprox(st.stress, EPSILON)
- && natoms()==st.natoms()
- && (std::fabs(eweight-st.eweight) < EPSILON)
- && (std::fabs(fweight-st.fweight) < EPSILON)
- && (std::fabs(sweight-st.sweight) < EPSILON)
- && (std::fabs(energy-st.energy) < EPSILON)
- ;
- if (!result) return result;
- for (size_t i=0;i<natoms();++i) {
- result = atoms[i]==st.atoms[i];
- if (!result) return result;
- }
- return result;
-}
-bool Structure::is_the_same(const Structure& st, double thr) const
-{
- bool result =
- cell.isApprox(st.cell, thr)
- && natoms()==st.natoms();
- if (!result) return result;
-
- size_t count=0;;
- for (size_t i=0;i<natoms();++i) {
- for (size_t j=i;j<natoms();++j) {
- result = atoms[i].is_the_same(st.atoms[j], thr);
- if (result) count++;
- }
- }
-
- return count==natoms() ? true : false;
-}
-int Structure::next_structure(std::ifstream &ifs) {
- std::string line;
- std::getline(ifs,line);
- if(line.empty()) return 0;
-
- std::getline(ifs,line);
- // the second line could be energy or
- // a scalling factors eweight fweight sweight
- std::stringstream stream(line);
- std::string temp;
- size_t count = 0;
- while(stream >> temp) { ++count;}
- // optional if second line is a weight
- if (count==3)
- std::getline(ifs,line);
-
- for (size_t i=0; i<6;++i)
- std::getline(ifs,line);
-
- int natoms=0;
- while(std::getline(ifs,line)) {
- if(line.empty()) break;
- if(line == "\r") break; // detects windows newline
- natoms++;
- }
- return natoms;
-}
-void Structure::clear_nn() {
- near_neigh_shift.clear();
- near_neigh_idx.clear();
-}
-
-std::vector<Atom>::iterator Structure::begin() {
- return atoms.begin();
-}
-
-std::vector<Atom>::iterator Structure::end() {
- return atoms.end();
-}
-
-std::vector<Atom>::const_iterator Structure::begin() const {
- return atoms.cbegin();
-}
-
-std::vector<Atom>::const_iterator Structure::end() const {
- return atoms.cend();
-}
-void Structure::dump_to_file(std::ostream& file, size_t prec) const {
- const int n = 5;
- file << label << std::endl;
- file << std::fixed << std::setprecision(prec);
- file << eweight << " " << fweight << " " << sweight << std::endl;
- file << energy << " " << T << std::endl;
-
- file
- << std::setw(prec+n) << cell(0,0) << " "
- << std::setw(prec+n) << cell(0,1) << " "
- << std::setw(prec+n) << cell(0,2) << " " << std::endl
- << std::setw(prec+n) << cell(1,0) << " "
- << std::setw(prec+n) << cell(1,1) << " "
- << std::setw(prec+n) << cell(1,2) << " " << std::endl
- << std::setw(prec+n) << cell(2,0) << " "
- << std::setw(prec+n) << cell(2,1) << " "
- << std::setw(prec+n) << cell(2,2) << " " << std::endl;
-
- file
- << std::setw(prec+n) << stress(0,0) << " "
- << std::setw(prec+n) << stress(0,1) << " "
- << std::setw(prec+n) << stress(0,2) << " " << std::endl
- << std::setw(prec+n) << stress(1,0) << " "
- << std::setw(prec+n) << stress(1,1) << " "
- << std::setw(prec+n) << stress(1,2) << " " << std::endl
- << std::setw(prec+n) << stress(2,0) << " "
- << std::setw(prec+n) << stress(2,1) << " "
- << std::setw(prec+n) << stress(2,2) << " " << std::endl;
-
- for (const auto& a : atoms) {
- file << std::setw(2) << a.symbol << " "
- << std::setw(prec+n) << a.position[0] << " "
- << std::setw(prec+n) << a.position[1] << " "
- << std::setw(prec+n) << a.position[2] << " "
- << std::setw(prec+n) << a.force[0] << " "
- << std::setw(prec+n) << a.force[1] << " "
- << std::setw(prec+n) << a.force[2] << std::endl;
- }
- file << std::endl;
-}
-void Structure::dump_to_file(const std::string& filepath, size_t prec) const {
- std::ofstream file(filepath, std::ios::app); // Open in append mode
- if (!file.is_open()) {
- std::cerr << "Error: Could not open file for writing: " << filepath << std::endl;
- return;
- }
- dump_to_file(file,prec);
- file.close();
-}
-std::set<Element> Structure::get_unique_elements() const{
- std::set<Element> unique_elements;
- for (const auto& a:atoms) unique_elements.insert(a);
- return unique_elements;
-}
diff --git a/src/structure_db.cpp b/src/structure_db.cpp
deleted file mode 100644
index f1aaa36f17971a77b3c18ab8a31e0653aab409fb..0000000000000000000000000000000000000000
--- a/src/structure_db.cpp
+++ /dev/null
@@ -1,356 +0,0 @@
-#include "tadah/core/element.h"
-#include <tadah/mlip/structure_db.h>
-#include <tadah/core/periodic_table.h>
-#include <cstdio>
-#include <cctype> // For std::isspace
-
-StructureDB::StructureDB() {
- PeriodicTable::initialize();
-}
-StructureDB::StructureDB(Config &config) {
- PeriodicTable::initialize();
- add(config);
-}
-
-void StructureDB::add(const std::string fn) {
- std::ifstream ifs(fn);
- if (!ifs.is_open()) {
- throw std::runtime_error("DBFILE does not exist: "+fn);
- }
- parseFile(fn);
- while (true) {
- structures.push_back(Structure());
- int t = structures.back().read(ifs);
-
- // did we read structure succesfully?
- // if not remove last object from the list
- if (t==1) structures.pop_back();
-
- if (ifs.eof()) break;
- }
- ifs.close();
-}
-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);
- }
-
- // move iterator to the first structure we want to read
- for (size_t i=0; i<first; ++i) {
- Structure::next_structure(ifs);
- }
-
- int count=0;
- while (count++<N) {
- structures.push_back(Structure());
- int t = structures.back().read(ifs);
-
- // did we read structure succesfully?
- // if not remove last object from the list
- if (t==1) structures.pop_back();
-
- if (ifs.eof()) break;
- }
- ifs.close();
-
- return --count;
-}
-
-void StructureDB::add(const Structure &s) {
- structures.push_back(s);
-}
-void StructureDB::add(const StructureDB &stdb) {
- for (const auto &s: stdb) add(s);
-}
-
-void StructureDB::remove(size_t i) {
- structures.erase(structures.begin()+i);
-}
-
-void StructureDB::add(Config &config) {
- for (const std::string &s : config("DBFILE")) {
- dbidx.push_back(size());
- add(s);
- }
- dbidx.push_back(size());
-}
-
-size_t StructureDB::size() const {
- return structures.size();
-}
-
-size_t StructureDB::size(size_t n) const {
- return dbidx[n+1]-dbidx[n];
-}
-
-Structure &StructureDB::operator()(size_t s) {
- return structures[s];
-}
-
-const Structure &StructureDB::operator()(size_t s) const {
- return structures[s];
-}
-
-Atom &StructureDB::operator()(size_t s, size_t a) {
- return structures[s].atoms[a];
-}
-size_t StructureDB::calc_natoms() const {
- size_t natoms=0;
- for (auto struc: structures) natoms += struc.natoms();
- return natoms;
-}
-size_t StructureDB::calc_natoms(size_t n) const {
- size_t start = dbidx[n];
- size_t stop = dbidx[n+1];
- size_t natoms=0;
- for (size_t i=start; i<stop; ++i) {
- natoms += (*this)(i).natoms();
- }
- return natoms;
-}
-std::set<Element> StructureDB::get_unique_elements() const {
- std::set<Element> s;
- for (const auto & st: structures) {
- std::set<Element> u = st.get_unique_elements();
- s.insert(u.cbegin(),u.cend());
- }
- return s;
-}
-std::set<Element> StructureDB::find_unique_elements(const std::string &fn) {
- std::set<Element> s;
- std::ifstream ifs(fn);
- std::string line;
-
- if (!ifs.is_open()) {
- std::cerr << "Could not open the file." << std::endl;
- }
-
- char symbol[3];
- while (std::getline(ifs, line)) {
- // the second line could be energy or
- // a scalling factors eweight fweight sweight
- std::getline(ifs,line);
- std::stringstream stream(line);
- size_t count = std::distance(std::istream_iterator<std::string>(stream),
- std::istream_iterator<std::string>());
-
- if (count == 3)
- std::getline(ifs,line);
-
- for (size_t i=0; i<6; ++i)
- std::getline(ifs,line);
-
- while (std::getline(ifs, line)) {
- if(line.empty()) break;
- sscanf(line.c_str(), "%2s", symbol);
- s.insert(PeriodicTable::find_by_symbol(symbol));
- }
-
- }
- return s;
-}
-std::set<Element> StructureDB::find_unique_elements(const Config &c) {
- std::set<Element> s;
- for (const auto& fn: c("DBFILE")) {
- std::set<Element> temp = find_unique_elements(fn);
- s.insert(temp.begin(), temp.end());
- }
- return s;
-}
-
-template <typename T,typename U>
-std::pair<T,U> operator+(const std::pair<T,U> &l,const std::pair<T,U> &r) {
- return {l.first+r.first,l.second+r.second};
-}
-std::pair<int,int> StructureDB::count(const Config &config) {
- std::pair<int, int> res=std::make_pair(0,0);
- for (const std::string &fn : config("DBFILE")) {
- res = res + count(fn);
- }
- return res;
-}
-
-std::pair<int,int> StructureDB::count(const std::string fn){
- int nstruc=0;
- int natoms_tot=0;
-
- std::ifstream ifs(fn);
- if (!ifs.is_open()) {
- throw std::runtime_error("DBFILE does not exist: "+fn);
- }
- int natoms = Structure::next_structure(ifs);
- while(natoms) {
- natoms_tot += natoms;
- natoms = Structure::next_structure(ifs);
- nstruc++;
- }
-
- std::pair<int,int> res = std::make_pair(nstruc,natoms_tot);
- return res;
-}
-void StructureDB::clear_nn() {
- for (auto &struc: structures) struc.clear_nn();
-}
-void StructureDB::check_atoms_key(Config &config, std::set<Element> &unique_elements) {
- bool error=false;
-
- if (config.exist("ATOMS")) {
- // user set this key so here we check does it correspond to unique_elements
- if (unique_elements.size()!=config.size("ATOMS")) {
- error=true;
- }
-
- auto set_it = unique_elements.begin();
- auto atoms_it = config("ATOMS").begin();
- while (set_it != unique_elements.end() && atoms_it != config("ATOMS").end()) {
- if (set_it->symbol != (*atoms_it)) {
- error=true;
- }
- ++set_it;
- ++atoms_it;
- }
- if (error) {
- throw std::runtime_error("\n"
- "Mismatch between elements in datasets and ATOMS in the config file.\n"
- "Please either update the ATOMS in the config file or remove ATOMS\n"
- "key completely. Tadah! will automatically configure this key.\n"
- );
- }
- } else {
- for (const auto &s : unique_elements) config.add("ATOMS", s.symbol);
- }
-}
-void StructureDB::check_watoms_key(Config &config, std::set<Element> &unique_elements) {
- bool error=false;
-
- if (config.exist("WATOMS")) {
- // user set this key so here we check does it correspond to unique_elements
- if (unique_elements.size()!=config.size("WATOMS")) {
- error=true;
- }
- if (error) {
- throw std::runtime_error("\n"
- "Mismatch between elements in datasets and WATOMS in the config file.\n"
- "Please either update the WATOMS in the config file or remove WATOMS\n"
- "key completely. In the latter case Tadah! will use default values.\n"
- );
- }
- } else {
- for (const auto &s : unique_elements) config.add("WATOMS", s.Z);
- }
-}
-
-std::vector<Structure>::iterator StructureDB::begin() {
- return structures.begin();
-}
-
-std::vector<Structure>::iterator StructureDB::end() {
- return structures.end();
-}
-
-std::vector<Structure>::const_iterator StructureDB::begin() const {
- return structures.cbegin();
-}
-
-std::vector<Structure>::const_iterator StructureDB::end() const {
- return structures.cend();
-}
-void StructureDB::dump_to_file(const std::string& filepath, size_t prec) const {
- std::ofstream file(filepath, std::ios::app); // Open in append mode
- if (!file.is_open()) {
- std::cerr << "Error: Could not open file for writing: " << filepath << std::endl;
- return;
- }
- for (const auto &s : structures) {
- s.dump_to_file(file,prec);
- }
- file.close();
-}
-
-std::string StructureDB::summary() const {
- std::string str = "# of structures : " + to_string(structures.size());
-
- str += " | # of atoms : " + to_string(calc_natoms());
-
- str += " | Elements : ";
- std::set<Element> ue = get_unique_elements();
- for (const auto &e: ue) str+= e.symbol + " ";
- str+="\n";
- return str;
-}
-void StructureDB::parseFile(const std::string& filename)
-{
- std::ifstream fin(filename, std::ios::in | std::ios::binary);
- if (!fin.is_open())
- {
- std::cerr << "Error: could not open file " << filename << "\n";
- return;
- }
-
- std::size_t header_size = 9;
-
- // Increase buffer size to speed up I/O on large files.
- static const size_t BUFSIZE = 100ULL << 20; // 100 MiB
- char* buffer = new char[BUFSIZE];
- fin.rdbuf()->pubsetbuf(buffer, BUFSIZE);
-
- std::vector<size_t> blockLineCounts;
- blockLineCounts.reserve(10000); // Pre-allocate to reduce repeated allocations
-
- size_t currentBlockCount = 0;
- std::string line;
-
- while (true)
- {
- if (!std::getline(fin, line))
- {
- // End of file or read error
- break;
- }
-
- if (isBlankLine(line))
- {
- // We reached the end of the current block
- if (currentBlockCount > 0)
- {
- blockLineCounts.push_back(currentBlockCount-header_size);
- currentBlockCount = 0;
- }
- }
- else
- {
- // Non-empty line => belongs to the current block
- currentBlockCount++;
- }
- }
-
- // If the last block didn’t end with a blank line, close it out
- if (currentBlockCount > 0)
- {
- blockLineCounts.push_back(currentBlockCount-header_size);
- }
-
- fin.close();
- delete[] buffer;
-
- // Print the results
- std::cout << "Found " << blockLineCounts.size() << " blocks.\n";
- for (size_t i = 0; i < blockLineCounts.size(); i+=1000)
- {
- std::cout << "Block " << i << " has "
- << blockLineCounts[i] << " atoms\n";
- }
-}
-
-bool StructureDB::isBlankLine(const std::string& line) const
-{
- for (char c : line)
- {
- if (!std::isspace(static_cast<unsigned char>(c)))
- {
- return false;
- }
- }
- return true;
-}