diff --git a/include/tadah/mlip/descriptors_calc.h b/include/tadah/mlip/descriptors_calc.h
index d4e9a77c918be467f1e873b1377aeb4d10996c11..3bb374319c8cec040f2241c64de9f5a5de911381 100644
--- a/include/tadah/mlip/descriptors_calc.h
+++ b/include/tadah/mlip/descriptors_calc.h
@@ -1,11 +1,13 @@
#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>
#include <tadah/mlip/st_descriptors.h>
#include <tadah/mlip/st_descriptors_db.h>
+#include <tadah/mlip/structure_neighbor_view.h>
#include <tadah/core/config.h>
#include <tadah/models/descriptors/d2/d2_base.h>
#include <tadah/models/descriptors/d3/d3_base.h>
@@ -108,7 +110,7 @@ class DescriptorsCalc: public DC_Base {
StDescriptorsDB calc(const StructureDB &st_db);
/** Calculate density vector for the structure */
- void calc_rho(const Structure &st, StDescriptors &std);
+ void calc_rho(const Structure &st, StructureNeighborView &st_nb, StDescriptors &std);
C2 c2;
C3 c3;
@@ -120,7 +122,7 @@ class DescriptorsCalc: public DC_Base {
template <typename T1, typename T2, typename T3>
DescriptorsCalc(Config &c,T1 &t1,T2 &t2,T3 &t3);
- void calc(const Structure &st, StDescriptors &std);
+ void calc(const Structure &st, 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 f7e876f565aa9a68dfdf80a712b0c052b68b7378..1c47da2a3d599856a2563c627ce095506ad9a1ec 100644
--- a/include/tadah/mlip/descriptors_calc.hpp
+++ b/include/tadah/mlip/descriptors_calc.hpp
@@ -4,6 +4,8 @@
#include <cstddef>
#include <tadah/mlip/descriptors_calc.h>
#include <tadah/core/periodic_table.h>
+#include <tadah/mlip/atom.h>
+#include <tadah/mlip/structure_neighbor_view.h>
#include <cstdio>
@@ -160,7 +162,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::common_constructor() {
config.add("DSIZE",dsize);
}
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, StDescriptors &st_d) {
+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;
size_t s = dm.rhoi_size()+dm.rhoip_size();
@@ -168,22 +170,38 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_rho(const Structure &st, StDescrip
rhos.set_zero();
Vec3d delij;
+ Vec3d rj;
for (size_t i=0; i<st.natoms(); ++i) {
const Atom &a1 = st(i);
- int Zi = a1.Z;
- for (size_t jj=0; jj<st.nn_size(i); ++jj) {
- const Vec3d &a2pos = st.nn_pos(i,jj);
+ 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);
+
+ size_t gj = NbPtr[jj];
+ size_t j = st.globalToLocal(gj); // j is local to st
+
+ auto sx = *st_nb.getShiftXPtr(j);
+ auto sy = *st_nb.getShiftYPtr(j);
+ auto sz = *st_nb.getShiftZPtr(j);
+
+ rj(0) = st.x(j) + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ 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);
+
//delij = a1.position - a2pos;
- delij[0] = a1.position[0] - a2pos[0];
- delij[1] = a1.position[1] - a2pos[1];
- delij[2] = a1.position[2] - a2pos[2];
+ delij[0] = a1.x() - rj(0);
+ delij[1] = a1.y() - rj(1);
+ delij[2] = a1.z() - rj(2);
//double rij_sq = delij * delij;
double rij_sq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
if (rij_sq > rcut_mb_sq) continue;
- size_t neighIdx = st.near_neigh_idx[i][jj];
- int Zj = st(neighIdx).Z;
+ // size_t neighIdx = st.near_neigh_idx[i][jj];
+ size_t neighIdx = st_nb.getAAIndex(i, j, jj);
+ int Zj = st(neighIdx).atomicNumber();
double rij = sqrt(rij_sq);
dm.calc_rho(Zi,Zj,rij,rij_sq,delij,st_d.get_rho(i));
}
@@ -191,7 +209,7 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc_rho(const Structure &st, StDescrip
}
template <typename D2, typename D3, typename DM, typename C2, typename C3, typename CM>
-void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, StDescriptors &st_d) {
+void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, const StructureNeighborView &st_nb, StDescriptors &st_d) {
bool init2b = config.get<bool>("INIT2B");
bool initmb = config.get<bool>("INITMB");
@@ -241,23 +259,38 @@ void DescriptorsCalc<D2,D3,DM,C2,C3,CM>::calc(const Structure &st, StDescriptors
bool use_stress = config.get<bool>("STRESS");
Vec3d delij;
+ Vec3d rj;
for (size_t i=0; i<st.natoms(); ++i) {
const Atom &a1 = st(i);
aed_type &aed = st_d.get_aed(i);
- int Zi = a1.Z;
- for (size_t jj=0; jj<st.nn_size(i); ++jj) {
- const Vec3d &a2pos = st.nn_pos(i,jj);
+ 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);
+
+ size_t gj = NbPtr[jj];
+ size_t j = st.globalToLocal(gj); // j is local to st
+
+ auto sx = *st_nb.getShiftXPtr(j);
+ auto sy = *st_nb.getShiftYPtr(j);
+ auto sz = *st_nb.getShiftZPtr(j);
+
+ rj(0) = st.x(j) + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ 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);
- delij[0] = a1.position[0] - a2pos[0];
- delij[1] = a1.position[1] - a2pos[1];
- delij[2] = a1.position[2] - a2pos[2];
+ delij[0] = a1.x() - rj(0);
+ delij[1] = a1.y() - rj(1);
+ delij[2] = a1.z() - rj(2);
double rij_sq = delij[0]*delij[0] + delij[1]*delij[1] + delij[2]*delij[2];
if (rij_sq > rcut_max_sq) continue;
- size_t neighIdx = st.near_neigh_idx[i][jj];
- int Zj = st(neighIdx).Z;
+ // size_t neighIdx = st.near_neigh_idx[i][jj];
+ size_t neighIdx = st_nb.getAAIndex(i, j, jj);
+ int Zj = st(neighIdx).atomicNumber();
double rij = sqrt(rij_sq);
double rij_inv = 1.0/rij;
diff --git a/include/tadah/mlip/design_matrix/design_matrix.h b/include/tadah/mlip/design_matrix/design_matrix.h
index a19ba9ba007267d820b593d084b768fa9f3b45d7..c7dcaa79627a007d3ac197e0eadafc894c413f9d 100644
--- a/include/tadah/mlip/design_matrix/design_matrix.h
+++ b/include/tadah/mlip/design_matrix/design_matrix.h
@@ -246,11 +246,11 @@ public:
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) {
ws->Tlabels(j) = 1;
- ws->T(j++) = a.force(0)*fscale;
+ ws->T(j++) = a.fx()*fscale;
ws->Tlabels(j) = 1;
- ws->T(j++) = a.force(1)*fscale;
+ ws->T(j++) = a.fy()*fscale;
ws->Tlabels(j) = 1;
- ws->T(j++) = a.force(2)*fscale;
+ ws->T(j++) = a.fz()*fscale;
}
}
if (stress) {
diff --git a/include/tadah/mlip/models/m_tadah_base.h b/include/tadah/mlip/models/m_tadah_base.h
index 4522eefdd94a1cc54ffe39b68cd3adb6a5699cf7..7ac57c8c5a46b0eed89810578baecd506a42b7e2 100644
--- a/include/tadah/mlip/models/m_tadah_base.h
+++ b/include/tadah/mlip/models/m_tadah_base.h
@@ -5,6 +5,7 @@
#include <tadah/mlip/structure_db.h>
#include <tadah/mlip/st_descriptors.h>
#include <tadah/mlip/st_descriptors_db.h>
+#include <tadah/mlip/structure_neighbor_view.h>
#include <tadah/mlip/normaliser.h>
#include <tadah/mlip/descriptors_calc_base.h>
#include <tadah/core/core_types.h>
@@ -40,15 +41,15 @@ public:
/** \brief Predict total force on an atom a. */
virtual void fpredict(const size_t a, force_type &v,
- const StDescriptors &std, const Structure &st);
+ 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,
- const Structure &st);
+ const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict energy, forces and stresses for the Structure st. */
virtual Structure predict(const Config &c,
- /*not const*/ StDescriptors &std, const Structure &st);
+ /*not const*/ StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
/** \brief Predict energy, forces and stresses for a set of Structures.
*
@@ -56,24 +57,24 @@ public:
*/
virtual StructureDB predict(const Config &c,
/*not const*/ StDescriptorsDB &st_desc_db,
- const StructureDB &stdb);
+ 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, DC_Base &dc);
+ virtual StructureDB predict(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);
+ virtual 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,
- const Structure &st);
+ 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);
+ 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,
- const StDescriptors &std, const Structure &st);
+ const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb);
/** Return potential file. */
@@ -98,10 +99,10 @@ public:
*/
virtual void train(StDescriptorsDB &, const StructureDB &) {};
- virtual StructureDB predict(Config config_pred, StructureDB &stdb, DC_Base &dc,
+ virtual StructureDB predict(Config config_pred, StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc,
aed_type &predicted_error)=0;
- virtual StructureDB predict(StructureDB &stdb)=0;
+ virtual StructureDB predict(StructureDB &stdb, const NeighborListDB &nldb)=0;
};
}
diff --git a/include/tadah/mlip/structure.h b/include/tadah/mlip/structure.h
index 3ebe47f34f8730cb83b38c8ce066d044d0f591e7..908f8180179ce5fc8822df5b85754d3e4321eb0c 100644
--- a/include/tadah/mlip/structure.h
+++ b/include/tadah/mlip/structure.h
@@ -7,6 +7,7 @@
#include <fstream>
#include <stdexcept>
#include <iterator> // for std::forward_iterator_tag
+#include <set>
// HPC types (Matrix3d, stress_type) come from here:
#include <tadah/core/periodic_table.h>
@@ -104,6 +105,14 @@ inline std::size_t globalToLocal(std::size_t gIdx) const {
return gIdx - offset_;
}
+ /** Return unique elements (atomic numbers) for all Structures. */
+ std::set<int> get_unique_elements() const {
+ std::set<int> unique_elements;
+ for (size_t a=0; a<natoms(); ++a) unique_elements.insert(atomicNumber(a));
+ return unique_elements;
+ }
+
+
// --------------------------------------------------
// Minimal I/O
// --------------------------------------------------
diff --git a/include/tadah/mlip/structure_db.h b/include/tadah/mlip/structure_db.h
index 07d2c23b99af4e4c84cb274bbfb9567a1e88aae2..89090bacd9f70b37326e3d353baea3734c68f672 100644
--- a/include/tadah/mlip/structure_db.h
+++ b/include/tadah/mlip/structure_db.h
@@ -7,6 +7,7 @@
#include <fstream>
#include <iostream>
#include <sstream>
+#include <set>
#include <cctype>
#include <stdexcept>
@@ -65,6 +66,16 @@ public:
std::vector<std::vector<std::size_t>> atomsPerStructurePerFile;
};
+
+ std::set<int> get_unique_elements() const {
+ std::set<int> s;
+ for (const auto & atomic_num : atomicNumber_) {
+ s.insert(atomic_num);
+ }
+ return s;
+ }
+
+
private:
/**
* @brief HEADERSIZE_ references the 9 lines forming a header in each structure file block.
diff --git a/include/tadah/mlip/structure_neighbor_view.h b/include/tadah/mlip/structure_neighbor_view.h
index fb871478d2b54c0ab262d3475db2210b5cd7d0b7..a8ede89daf410b64a933146742c851dcee3c3051 100644
--- a/include/tadah/mlip/structure_neighbor_view.h
+++ b/include/tadah/mlip/structure_neighbor_view.h
@@ -123,17 +123,6 @@ public:
return globalList_.getShiftZPtr(gIdx);
}
-private:
- /**
- * @brief Reference to the single global HPC neighbor list (covering all atoms).
- */
- const NeighborList &globalList_;
-
- /**
- * @brief The targeted structure from the HPC DB, which provides offset() and natoms().
- */
- const Structure &str_;
-
/**
* @brief Returns the local index (ii) of the i-th atom within the j-th atom's neighbor list.
*
@@ -165,7 +154,7 @@ private:
* @param jj Position of j in i's neighbor list.
* @return The local index (ii) of i in j's neighbor list.
*/
- inline size_t getAAIndex(const size_t i, const size_t j, const size_t jj)
+ inline size_t getAAIndex(const size_t i, const size_t j, const size_t jj) const
{
// Verify bounds on i, j, and jj to avoid out-of-range errors.
// i, j, and jj must be valid for this structure.
@@ -199,6 +188,18 @@ private:
return ii;
}
+
+private:
+ /**
+ * @brief Reference to the single global HPC neighbor list (covering all atoms).
+ */
+ const NeighborList &globalList_;
+
+ /**
+ * @brief The targeted structure from the HPC DB, which provides offset() and natoms().
+ */
+ const Structure &str_;
+
};
} // namespace mlip
diff --git a/src/dm_bf_linear.cpp b/src/dm_bf_linear.cpp
index 77f84db27b2015bdc5aed340ec33ebab4c2a7f7c..43de0a4da9cb1ca22030eef0656ff41be8b03e98 100644
--- a/src/dm_bf_linear.cpp
+++ b/src/dm_bf_linear.cpp
@@ -26,15 +26,14 @@ void DM_BF_Linear::calc_phi_force_rows(phi_type &Phi, size_t &row,
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);
+ 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);
+ 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)
+ 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];
for (size_t k = 0; k < 3; ++k) {
@@ -52,24 +51,36 @@ void DM_BF_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
const StructureNeighborView &st_nb,
const StDescriptors &st_d) {
double V_inv = 1 / st.get_volume();
+ Vec3d rj;
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 Vec3dSoAConstView ri = st.positionView(a);
- const aed_type& aedi = st_d.get_aed(a);
- for (size_t jj=0; jj<NNbrs; ++jj) {
+ Vec3dSoAConstView ra = 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];
+ size_t j = st.globalToLocal(gj); // j is local to st
+ size_t aa = st_nb.getAAIndex(a, j, jj);
+
+ auto sx = *st_nb.getShiftXPtr(j);
+ auto sy = *st_nb.getShiftYPtr(j);
+ auto sz = *st_nb.getShiftZPtr(j);
+
+ rj(0) = st.x(j) + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ 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 Vec3d &rj = st.nn_pos(a, jj);
+ const fd_type &fdij = st_d.fd[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) {
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));
+ fac[mn] * (ra(x) - rj(x));
}
mn++;
}
diff --git a/src/dm_bf_polynomial2.cpp b/src/dm_bf_polynomial2.cpp
index b139f79e51b172244ae9e8d67feff18a8b468c2c..cf7dbd5cb955359a58f56023f6944ee1b5720088 100644
--- a/src/dm_bf_polynomial2.cpp
+++ b/src/dm_bf_polynomial2.cpp
@@ -40,21 +40,17 @@ void DM_BF_Polynomial2::calc_phi_force_rows(phi_type &Phi,
{
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 ra = st.positionView(a);
+ Vec3dSoAConstView ra = 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 aed_type& aedi = st_d.get_aed(a);
+ 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);
-
for (size_t k=0; k<3; ++k) {
size_t b=0;
for (size_t i=0; i<fdij.rows(); ++i) {
@@ -77,25 +73,25 @@ void DM_BF_Polynomial2::calc_phi_stress_rows(phi_type &Phi,
const StDescriptors &st_d)
{
double V_inv = 1 / st.get_volume();
+ 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);
- const Vec3dSoAView ra = st.positionView(a);
+ Vec3dSoAConstView 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);
+ size_t aa = st_nb.getAAIndex(a, j, jj);
- auto sx = st.getShiftXPtr(j);
- auto sy = st.getShiftYPtr(j);
- auto sz = st.getShiftZPtr(j);
+ auto sx = *st_nb.getShiftXPtr(j);
+ auto sy = *st_nb.getShiftYPtr(j);
+ auto sz = *st_nb.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);
+ rj(0) = st.x(j) + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ 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];
diff --git a/src/dm_kern_base.cpp b/src/dm_kern_base.cpp
index 28df869f4760ae21bc912834fe8c904cc58fa48f..85d62718e92237a2a2a598824bafb5e703547c42 100644
--- a/src/dm_kern_base.cpp
+++ b/src/dm_kern_base.cpp
@@ -28,13 +28,17 @@ void DM_Kern_Base::calc_phi_energy_row(phi_type &Phi, size_t &row, const double
}
void DM_Kern_Base::calc_phi_force_rows(phi_type &Phi, size_t &row, const double fac,
const Structure &st, const StructureNeighborView &st_nb, const StDescriptors &st_d) {
- for (size_t i=0; i<st.natoms(); ++i) {
- const aed_type& aedi = st_d.get_aed(i);
- for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
- size_t j=st.near_neigh_idx[i][jj];
- size_t ii = st.get_nn_iindex(i,j,jj);
- const fd_type &fdji = st_d.fd[j][ii];
- const fd_type &fdij = st_d.fd[i][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);
+ Vec3dSoAConstView ra = 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 = 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);
for (size_t b=0; b<basis.cols(); ++b) {
for (size_t k=0; k<3; ++k) {
@@ -50,21 +54,34 @@ void DM_Kern_Base::calc_phi_stress_rows(phi_type &Phi, size_t &row, const doubl
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;
- const aed_type& aedi = st_d.get_aed(i);
- for (size_t jj=0; jj<st_d.fd[i].size(); ++jj) {
- size_t j=st.near_neigh_idx[i][jj];
- size_t ii = st.get_nn_iindex(i,j,jj);
- const fd_type &fdji = st_d.fd[j][ii];
- const fd_type &fdij = st_d.fd[i][jj];
+ 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);
+ 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);
+
+ auto sx = *st_nb.getShiftXPtr(j);
+ auto sy = *st_nb.getShiftYPtr(j);
+ auto sz = *st_nb.getShiftZPtr(j);
+
+ rj(0) = st.x(j) + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ 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 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 b=0; b<basis.cols(); ++b) {
- Phi(row+mn,b) += V_inv*0.5*fac[mn]*(ri(x)-rj(x))*
+ Phi(row+mn,b) += V_inv*0.5*fac[mn]*(ra(x)-rj(x))*
((*this).prime(basis.col(b),aedi,fdij(y)) -
(*this).prime(basis.col(b),aedj,fdji(y)));
}
diff --git a/src/dm_kern_linear.cpp b/src/dm_kern_linear.cpp
index dc5a89836ffc7bea4d64a7fc3ecd1cc37aec1b6a..2a675e5b9349074e28409d54953bc92541c3a3d9 100644
--- a/src/dm_kern_linear.cpp
+++ b/src/dm_kern_linear.cpp
@@ -32,9 +32,12 @@ void DM_Kern_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)
{
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);
+ const std::size_t* NbPtr = st_nb.getNeighborsPtr(a);
+ size_t NNbrs = st_nb.numNeighbors(a);
+ for (size_t jj=0; jj<NNbrs; ++jj) {
+ const size_t gj = NbPtr[jj];
+ 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)-
st_d.fd[j][aa](k))*fac;
@@ -51,18 +54,30 @@ void DM_Kern_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);
+ 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);
+ 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];
+
+ auto sx = *st_nb.getShiftXPtr(j);
+ auto sy = *st_nb.getShiftYPtr(j);
+ auto sz = *st_nb.getShiftZPtr(j);
+
+ rj(0) = st.x(j) + sx*st.cell(0,0) + sy*st.cell(0,1) + sz*st.cell(0,2);
+ 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);
+
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]*(ri(x)-rj(x));
+ 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];
}
@@ -73,5 +88,6 @@ void DM_Kern_Linear::calc_phi_stress_rows(phi_type &Phi, size_t &row,
}
row += 6;
}
+
}
}
diff --git a/src/lammps_structure_writer.cpp b/src/lammps_structure_writer.cpp
index 87f14f57cc91aafb73f2a5207bf1ae92ec233511..86f5468acc3bc30395a9e0a70a73580999ce3f7d 100644
--- a/src/lammps_structure_writer.cpp
+++ b/src/lammps_structure_writer.cpp
@@ -21,19 +21,21 @@ void LammpsStructureWriter::write_data(const StructureDB &stdb, const std::strin
const Structure &st = stdb(structIndex);
// compute number of atoms for a given element
- const auto &elements = st.get_unique_elements();
+ const auto &atomicNumbers = st.get_unique_elements();
std::map<std::string, size_t> nelements;
std::map<std::string, size_t> type;
// Initialize the count for each element
- for (const auto& element : elements) {
- nelements[element.symbol] = 0;
+ for (const auto& atomicNum : atomicNumbers) {
+ std::string symbol = PeriodicTable::find_by_Z(atomicNum).symbol;
+ nelements[symbol] = 0;
}
// then count...
size_t t=0;
for (const auto &atom : st) {
- nelements[atom.symbol]++;
- if (type.find(atom.symbol) == type.end()) {
- type[atom.symbol] = ++t;
+ std::string symbol = PeriodicTable::find_by_Z(atom.atomicNumber()).symbol;
+ nelements[symbol]++;
+ if (type.find(symbol) == type.end()) {
+ type[symbol] = ++t;
}
}
@@ -74,14 +76,15 @@ 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;
file << std::right << std::fixed << std::setw(w) << idx;
- file << std::right << std::fixed << std::setw(w) << type[atom.symbol];
+ file << std::right << std::fixed << std::setw(w) << type[symbol];
file << std::right << std::fixed << std::setw(w)
- << std::setprecision(p) << atom.position(0);
+ << std::setprecision(p) << atom.x();
file << std::right << std::fixed << std::setw(w)
- << std::setprecision(p) << atom.position(1);
+ << std::setprecision(p) << atom.y();
file << std::right << std::fixed << std::setw(w)
- << std::setprecision(p) << atom.position(2);
+ << std::setprecision(p) << atom.z();
file << std::endl;
idx++;
}
diff --git a/src/m_tadah_base.cpp b/src/m_tadah_base.cpp
index b1ea7af6adfd20d81a208879d6c06566ca0ffe21..51fa45389ac04ec74a7988eabd95b073a982bff7 100644
--- a/src/m_tadah_base.cpp
+++ b/src/m_tadah_base.cpp
@@ -4,7 +4,7 @@ namespace tadah {
namespace mlip {
void M_Tadah_Base::
fpredict(const size_t a, force_type &v,
- const StDescriptors &std, const Structure &st)
+ const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
for (size_t jj=0; jj<st.nn_size(a); ++jj) {
size_t j=st.near_neigh_idx[a][jj];
@@ -18,10 +18,10 @@ fpredict(const size_t a, force_type &v,
}
}
force_type M_Tadah_Base::
-fpredict(const size_t a, const StDescriptors &std, const Structure &st)
+fpredict(const size_t a, const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
force_type v;
- fpredict(a,v,std,st);
+ fpredict(a,v,std,st, st_nb);
return v;
}
double M_Tadah_Base::
@@ -35,7 +35,7 @@ epredict(const StDescriptors &std)
}
Structure M_Tadah_Base::
-predict(const Config &c, StDescriptors &std, const Structure &st)
+predict(const 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)
@@ -45,23 +45,23 @@ predict(const Config &c, StDescriptors &std, const Structure &st)
if (c.get<bool>("FORCE") && !c.get<bool>("STRESS")) {
for (size_t a=0; a<st.natoms(); ++a) {
- st_.atoms[a].force = fpredict(a,std,st);
+ st_.atoms[a].force = fpredict(a,std,st, st_nb);
}
}
else if (c.get<bool>("STRESS") && !c.get<bool>("FORCE")) {
- st_.stress = spredict(std,st);
+ st_.stress = spredict(std,st, st_nb);
}
else if (c.get<bool>("STRESS") && c.get<bool>("FORCE")) {
- stress_force_predict(std, st_);
+ stress_force_predict(std, st_, st_nb);
}
return st_;
}
StructureDB M_Tadah_Base::
-predict(const Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb)
+predict(const Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb, const NeighborListDB &nldb)
{
StructureDB stdb_;
stdb_.structures.resize(stdb.size());
@@ -69,13 +69,14 @@ predict(const Config &c, StDescriptorsDB &st_desc_db, const StructureDB &stdb)
#pragma omp parallel for
#endif
for (size_t i=0; i<stdb.size(); ++i) {
- stdb_(i) = predict(c,st_desc_db(i), stdb(i));
+ StructureNeighborView st_nb(nldb.nlist(),stdb(i));
+ stdb_(i) = predict(c,st_desc_db(i), stdb(i), st_nb);
}
return stdb_;
}
StructureDB M_Tadah_Base::
-predict(Config &c, const StructureDB &stdb, DC_Base &dc)
+predict(Config &c, const StructureDB &stdb, const NeighborListDB &nldb, DC_Base &dc)
{
StructureDB stdb_;
stdb_.structures.resize(stdb.size());
@@ -84,32 +85,33 @@ predict(Config &c, const StructureDB &stdb, DC_Base &dc)
#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));
+ StructureNeighborView st_nb(nldb.nlist(),stdb(i));
+ stdb_(i) = predict(c,st_d, stdb(i), st_nb);
}
return stdb_;
}
stress_type M_Tadah_Base::
-spredict(const StDescriptors &std, const Structure &st)
+spredict(const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
stress_type s;
s.set_zero();
for (size_t a=0; a<st.natoms(); ++a) {
- s += spredict(a,std,st);
+ s += spredict(a,std,st,st_nb);
}
return s;
}
stress_type M_Tadah_Base::
-spredict(const size_t a, const StDescriptors &std, const Structure &st)
+spredict(const size_t a, const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
stress_type s;
s.set_zero();
- spredict(a,s,std,st);
+ spredict(a,s,std,st,st_nb);
return s;
}
void M_Tadah_Base::
spredict(const size_t a, stress_type &s,
- const StDescriptors &std, const Structure &st)
+ const StDescriptors &std, const Structure &st, const StructureNeighborView &st_nb)
{
double V_inv = 1/st.get_volume();
const Vec3d &ri = st.atoms[a].position;
@@ -136,7 +138,7 @@ spredict(const size_t a, stress_type &s,
}
void M_Tadah_Base::
-stress_force_predict(const StDescriptors &std, Structure &st_)
+stress_force_predict(const StDescriptors &std, Structure &st_, const StructureNeighborView &st_nb)
{
stress_type s;
s.set_zero();