diff --git a/include/tadah/models/descriptors/dm/dm_base.h b/include/tadah/models/descriptors/dm/dm_base.h
index 276404efe6480afe62e7230a9415170a96e424b1..03ca07240c126db60eaee5062f157ec712518287 100644
--- a/include/tadah/models/descriptors/dm/dm_base.h
+++ b/include/tadah/models/descriptors/dm/dm_base.h
@@ -11,6 +11,7 @@
class DM_Base: public D_Base {
protected:
+ size_t rfidx=0;
size_t rhoisize=0;
public:
DM_Base() {};
@@ -125,5 +126,7 @@ public:
return s;
}
virtual std::string label()=0;
+ virtual size_t get_rfidx();
+ virtual void set_rfidx(size_t fidx_);
};
#endif
diff --git a/include/tadah/models/descriptors/dm/dm_mEAD.hpp b/include/tadah/models/descriptors/dm/dm_mEAD.hpp
index 2ff8f7282a1c4718759e16c447b145663d11e3ba..9c1dc124bcb10f3f626d3a0d5223b61109f067e0 100644
--- a/include/tadah/models/descriptors/dm/dm_mEAD.hpp
+++ b/include/tadah/models/descriptors/dm/dm_mEAD.hpp
@@ -56,7 +56,7 @@ void DM_mEAD<F>::calc_aed(
aed_type &aed)
{
size_t Lshift = 0;
- size_t ii=0;
+ size_t ii=rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const double f = fac[L][o];
@@ -92,7 +92,7 @@ void DM_mEAD<F>::calc_dXijdri_dXjidri(
double fc_ij = fcut->calc(rij);
double fcp_ij = fcut->calc_prime(rij);
- size_t ii=rhoisize;
+ size_t ii=rhoisize+rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const int lx = orbitals[L][o][0];
@@ -185,7 +185,7 @@ void DM_mEAD<F>::calc_dXijdri(
size_t Lshift = 0;
- size_t ii=rhoisize;
+ size_t ii=rhoisize+rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const int lx = orbitals[L][o][0];
@@ -258,7 +258,7 @@ void DM_mEAD<F>::calc_rho(
rho_type& rhoi,
const double scale)
{
- size_t ii=0;
+ size_t ii=rfidx;
double fc_ij = scale*weights[Zj]*fcut->calc(rij);
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
diff --git a/include/tadah/models/descriptors/dm/dm_mjoin.h b/include/tadah/models/descriptors/dm/dm_mjoin.h
index 5dd7656864d6c17cb10787adde030af29ee020cc..164adf93cae87d0a954dd8a73e4fafb0068aec9f 100644
--- a/include/tadah/models/descriptors/dm/dm_mjoin.h
+++ b/include/tadah/models/descriptors/dm/dm_mjoin.h
@@ -85,6 +85,7 @@ public:
rho_type& rho,
const double scale=1) override;
void set_fidx(size_t fidx_) override;
+ void set_rfidx(size_t rfidx_) override;
void init() override;
};
diff --git a/src/dm_base.cpp b/src/dm_base.cpp
index a76605da06946d74afac69f10e502df2037d6e3d..7fe63ac8afe5d8129e9ba4d86a7cdc63fff294a1 100644
--- a/src/dm_base.cpp
+++ b/src/dm_base.cpp
@@ -1 +1,3 @@
#include <tadah/models/descriptors/dm/dm_base.h>
+void DM_Base::set_rfidx(size_t rfidx_) { rfidx=rfidx_; }
+size_t DM_Base::get_rfidx() { return rfidx; }
diff --git a/src/dm_blip.cpp b/src/dm_blip.cpp
index 4748b6960530e80efb275052dadf208fa53965aa..f75dc44228e7660c5081f4ff686b132abae81b34 100644
--- a/src/dm_blip.cpp
+++ b/src/dm_blip.cpp
@@ -47,7 +47,7 @@ void DM_Blip::calc_aed(
aed_type &aed)
{
size_t Lshift = 0;
- size_t ii=0;
+ size_t ii=rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const double f = fac[L][o];
@@ -82,7 +82,7 @@ void DM_Blip::calc_dXijdri_dXjidri(
size_t Lshift = 0;
- size_t ii=rhoisize;
+ size_t ii=rhoisize+rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const int lx = orbitals[L][o][0];
@@ -173,7 +173,7 @@ void DM_Blip::calc_dXijdri(
size_t Lshift = 0;
- size_t ii=rhoisize;
+ size_t ii=rhoisize+rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const int lx = orbitals[L][o][0];
@@ -243,7 +243,7 @@ void DM_Blip::calc_rho(
rho_type& rhoi,
const double scale)
{
- size_t ii=0;
+ size_t ii=rfidx;
double fc_ij = scale*weights[Zj]*fcut->calc(rij);
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
diff --git a/src/dm_ead.cpp b/src/dm_ead.cpp
index 08bd5c5ee80443d5b2de52afb839d13278eec619..c241da9d918f5f2dda56296c65965c9452de3a11 100644
--- a/src/dm_ead.cpp
+++ b/src/dm_ead.cpp
@@ -54,7 +54,7 @@ void DM_EAD::calc_aed(
aed_type &aed)
{
size_t Lshift = 0;
- size_t ii=0;
+ size_t ii=rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const double f = fac[L][o];
@@ -89,7 +89,7 @@ void DM_EAD::calc_dXijdri_dXjidri(
size_t Lshift = 0;
- size_t ii=rhoisize;
+ size_t ii=rhoisize+rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const int lx = orbitals[L][o][0];
@@ -181,7 +181,7 @@ void DM_EAD::calc_dXijdri(
size_t Lshift = 0;
- size_t ii=rhoisize;
+ size_t ii=rhoisize+rfidx;
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
const int lx = orbitals[L][o][0];
@@ -251,7 +251,7 @@ void DM_EAD::calc_rho(
rho_type& rhoi,
const double scale)
{
- size_t ii=0;
+ size_t ii=rfidx;
double fc_ij = scale*weights[Zj]*fcut->calc(rij);
for (int L=0; L<=Lmax; ++L) {
for (size_t o=0; o<orbitals[L].size(); ++o) {
diff --git a/src/dm_eam.cpp b/src/dm_eam.cpp
index e8878e9f4583c82e5ace623a051cccb02cc0ecf5..94c0c48dfe3d901e81bb8a4932859e6fc26d0ab7 100644
--- a/src/dm_eam.cpp
+++ b/src/dm_eam.cpp
@@ -39,7 +39,7 @@ void DM_EAM::calc_aed(
aed_type &aed)
{
//std::cout << "rho[i]: " << rho(0) << std::endl;
- double p = rho(0)*ef.rdrho + 1.0;
+ double p = rho(rfidx)*ef.rdrho + 1.0;
int m = static_cast<int> (p);
m = std::max(1,std::min(m,ef.nrho-1));
p -= m;
@@ -47,10 +47,10 @@ void DM_EAM::calc_aed(
//std::vector<double> coeff = frho_spline[m];
double phi = ((frho_spline[m][3]*p + frho_spline[m][4])*p + frho_spline[m][5])*p + frho_spline[m][6];
- rho(1) = (frho_spline[m][0]*p + frho_spline[m][1])*p + frho_spline[m][2]; // lammps fp[i]
+ rho(1+rfidx) = (frho_spline[m][0]*p + frho_spline[m][1])*p + frho_spline[m][2]; // lammps fp[i]
//std::cout << "fp[i]: " << rho(1) << std::endl;
- if (rho(0) > ef.rhomax) phi += rho(1) * (rho(0)-ef.rhomax);
+ if (rho(rfidx) > ef.rhomax) phi += rho(1+rfidx) * (rho(rfidx)-ef.rhomax);
//phi *= scale[type[i]][type[i]];
aed(fidx) += phi;
@@ -80,7 +80,7 @@ void DM_EAM::calc_dXijdri_dXjidri(
// here rho_i[1] is a derivative of rho_i at local atom, similarly for rho_j
//double psip = rho_i[1]*rhojp + rho_j[1]*rhoip;
- double psip = rhoi(1)*rhojp + rhoj(1)*rhoip;
+ double psip = rhoi(rfidx+1)*rhojp + rhoj(rfidx+1)*rhoip;
double f = scale*psip*rinv;
fd_ij(fidx,0) = f*vec_ij[0];
@@ -107,7 +107,7 @@ void DM_EAM::calc_dXijdri(
double rhojp = (rhor_spline[m][0]*p + rhor_spline[m][1])*p + rhor_spline[m][2];
// here rho_i[1] is a derivative of rho_i at local atom
- double psip = rhoi(1)*rhojp;
+ double psip = rhoi(rfidx+1)*rhojp;
double f = scale*psip*rinv;
fd_ij(fidx,0) = f*vec_ij[0];
@@ -207,7 +207,7 @@ void DM_EAM::calc_rho(
p -= m;
p = std::min(p,1.0);
std::vector<double> coeff = rhor_spline[m];
- rho(0) += scale*( ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
+ rho(rfidx) += scale*( ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
}
void DM_EAM::init() {
diff --git a/src/dm_mjoin.cpp b/src/dm_mjoin.cpp
index c9e3eaea382171b8a98e5a837683834b3b50ad13..7e6a4023b922bedf9235b5061fa0afdd55ac40ec 100644
--- a/src/dm_mjoin.cpp
+++ b/src/dm_mjoin.cpp
@@ -18,15 +18,12 @@ DM_mJoin::DM_mJoin(Config &c) : DM_Base(c) {
s += ds.back()->size();
rhoisize += ds.back()->rhoi_size();
}
+ set_rfidx(0);
}
void DM_mJoin::calc_aed(rho_type& rho, aed_type &aed) {
- size_t rho_fidx = 0;
for (auto d : ds) {
- size_t rho_size = 2 * d->rhoi_size();
- rho_type rho_ptr(&rho[rho_fidx], rho_size, 1);
- d->calc_aed(rho_ptr, aed);
- rho_fidx += rho_size;
+ d->calc_aed(rho, aed);
}
}
void DM_mJoin::calc_dXijdri_dXjidri(
@@ -39,15 +36,10 @@ void DM_mJoin::calc_dXijdri_dXjidri(
rho_type& rhoj,
fd_type &fd_ij,
const double scale) {
- size_t rho_fidx = 0;
for (auto d : ds) {
- size_t rho_size = 2 * d->rhoi_size();
if (d->is_init_for_atoms(Zi,Zj) && rij < d->get_rcut()) {
- rho_type rhoi_ptr(&rhoi[rho_fidx],rho_size);
- rho_type rhoj_ptr(&rhoj[rho_fidx],rho_size);
- d->calc_dXijdri_dXjidri(Zi,Zj,rij,rij_sq,vec_ij,rhoi_ptr,rhoj_ptr,fd_ij,scale);
+ d->calc_dXijdri_dXjidri(Zi,Zj,rij,rij_sq,vec_ij,rhoi,rhoj,fd_ij,scale);
}
- rho_fidx += rho_size;
}
}
void DM_mJoin::calc_dXijdri(
@@ -59,14 +51,10 @@ void DM_mJoin::calc_dXijdri(
rho_type& rhoi,
fd_type &fd_ij,
const double scale) {
- size_t rho_fidx = 0;
for (auto d : ds) {
- size_t rho_size = 2 * d->rhoi_size();
if (d->is_init_for_atoms(Zi,Zj) && rij < d->get_rcut()) {
- rho_type rhoi_ptr(&rhoi[rho_fidx],rho_size);
- d->calc_dXijdri(Zi,Zj,rij,rij_sq,vec_ij,rhoi_ptr,fd_ij,scale);
+ d->calc_dXijdri(Zi,Zj,rij,rij_sq,vec_ij,rhoi,fd_ij,scale);
}
- rho_fidx += rho_size;
}
}
void DM_mJoin::init_rhoi(rho_type& rhoi) {
@@ -81,14 +69,10 @@ void DM_mJoin::calc_rho(
const Vec3d &vec_ij,
rho_type& rho,
const double scale) {
- size_t rho_fidx = 0;
for (auto d : ds) {
- size_t rho_size = 2 * d->rhoi_size();
if (d->is_init_for_atoms(Zi,Zj) && rij < d->get_rcut()) {
- rho_type rho_ptr(&rho[rho_fidx],rho_size);
- d->calc_rho(Zi,Zj,rij,rij_sq,vec_ij,rho_ptr,scale);
+ d->calc_rho(Zi,Zj,rij,rij_sq,vec_ij,rho,scale);
}
- rho_fidx += rho_size;
}
}
std::string DM_mJoin::label() {
@@ -96,6 +80,7 @@ std::string DM_mJoin::label() {
}
void DM_mJoin::set_fidx(size_t fidx_) {
+ fidx = fidx_;
ds[0]->set_fidx(fidx_);
size_t s = 0;
for (size_t i = 1; i < ds.size(); ++i) {
@@ -103,6 +88,16 @@ void DM_mJoin::set_fidx(size_t fidx_) {
ds[i]->set_fidx(fidx_ + s);
}
}
+void DM_mJoin::set_rfidx(size_t rfidx_) {
+ rfidx = rfidx_;
+ ds[0]->set_rfidx(rfidx_);
+ size_t s = 0;
+ for (size_t i = 1; i < ds.size(); ++i) {
+ s += ds[i - 1]->rhoi_size();
+ s += ds[i - 1]->rhoip_size();
+ ds[i]->set_rfidx(rfidx_ + s);
+ }
+}
void DM_mJoin::init() {
// not really sure about this