From ac898e635e300faef1ee243bed0695e639bcaac6 Mon Sep 17 00:00:00 2001
From: mkirsz <s1351949@sms.ed.ac.uk>
Date: Sun, 22 Dec 2024 07:34:49 +0000
Subject: [PATCH] each D-calculator to check(Zi,Zj) and cutoff

---
 include/tadah/models/descriptors/dm/dm_mEAD.hpp |  3 +++
 src/d2_blip.cpp                                 |  3 +++
 src/d2_bp.cpp                                   |  3 +++
 src/d2_eam.cpp                                  | 15 +++++++++------
 src/d2_lj.cpp                                   |  3 +++
 src/d2_mie.cpp                                  |  3 +++
 src/d2_mjoin.cpp                                |  6 ------
 src/d2_zbl.cpp                                  |  3 +++
 src/dm_blip.cpp                                 |  4 ++++
 src/dm_ead.cpp                                  |  3 +++
 src/dm_eam.cpp                                  | 15 +++++++++------
 src/dm_mjoin.cpp                                | 10 ++--------
 12 files changed, 45 insertions(+), 26 deletions(-)

diff --git a/include/tadah/models/descriptors/dm/dm_mEAD.hpp b/include/tadah/models/descriptors/dm/dm_mEAD.hpp
index 9c1dc12..c28018a 100644
--- a/include/tadah/models/descriptors/dm/dm_mEAD.hpp
+++ b/include/tadah/models/descriptors/dm/dm_mEAD.hpp
@@ -81,6 +81,7 @@ void DM_mEAD<F>::calc_dXijdri_dXjidri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   const double x = vec_ij[0];
   const double y = vec_ij[1];
   const double z = vec_ij[2];
@@ -175,6 +176,7 @@ void DM_mEAD<F>::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   const double x = vec_ij[0];
   const double y = vec_ij[1];
   const double z = vec_ij[2];
@@ -258,6 +260,7 @@ void DM_mEAD<F>::calc_rho(
   rho_type& rhoi,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t ii=rfidx;
   double fc_ij = scale*weights[Zj]*fcut->calc(rij);
   for (int L=0; L<=Lmax; ++L) {
diff --git a/src/d2_blip.cpp b/src/d2_blip.cpp
index 83bfc3d..7cf352f 100644
--- a/src/d2_blip.cpp
+++ b/src/d2_blip.cpp
@@ -51,6 +51,7 @@ void D2_Blip::calc_aed(
   aed_type &aed,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t i=fidx;
   double fc_ij = fcut->calc(rij);
   for (size_t c=0; c<mius.size(); c++) {
@@ -67,6 +68,7 @@ void D2_Blip::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t i=fidx;
   double fc_ij = fcut->calc(rij);
   double fcp_ij = fcut->calc_prime(rij);
@@ -83,6 +85,7 @@ void D2_Blip::calc_all(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
 
   size_t i=fidx;
   double fc_ij = fcut->calc(rij);
diff --git a/src/d2_bp.cpp b/src/d2_bp.cpp
index 4a980db..383995c 100644
--- a/src/d2_bp.cpp
+++ b/src/d2_bp.cpp
@@ -48,6 +48,7 @@ void D2_BP::calc_aed(
   aed_type &aed,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
 
   size_t i=fidx;
   double fc_ij = fcut->calc(rij);
@@ -64,6 +65,7 @@ void D2_BP::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t i=fidx;
   double fc_ij = fcut->calc(rij);
   double fcp_ij = fcut->calc_prime(rij);
@@ -80,6 +82,7 @@ void D2_BP::calc_all(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t i=fidx;
   double fc_ij = fcut->calc(rij);
   double fcp_ij = fcut->calc_prime(rij);
diff --git a/src/d2_eam.cpp b/src/d2_eam.cpp
index 18470f6..c168631 100644
--- a/src/d2_eam.cpp
+++ b/src/d2_eam.cpp
@@ -34,13 +34,14 @@ D2_EAM::D2_EAM(Config &c): D2_Base(c)
 }
 
 void D2_EAM::calc_aed(
-  const int ,
-  const int ,
+  const int Zi,
+  const int Zj,
   const double rij,
   const double,
   aed_type &aed,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r = rij;
   const double recip = 1.0/r;
   double p = r*ef.rdr + 1.0;
@@ -52,13 +53,14 @@ void D2_EAM::calc_aed(
   aed(fidx) += scale*z2*recip;
 }
 void D2_EAM::calc_dXijdri(
-  const int ,
-  const int ,
+  const int Zi,
+  const int Zj,
   const double rij,
   const double,
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   const double r = rij;
   const double recip = 1.0/r;
   double p = r*ef.rdr + 1.0;
@@ -76,14 +78,15 @@ void D2_EAM::calc_dXijdri(
 
 }
 void D2_EAM::calc_all(
-  const int ,
-  const int ,
+  const int Zi,
+  const int Zj,
   const double rij,
   const double,
   aed_type &aed,
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r = rij;
   const double recip = 1.0/r;
   double p = r*ef.rdr + 1.0;
diff --git a/src/d2_lj.cpp b/src/d2_lj.cpp
index b9369ef..64f33a3 100644
--- a/src/d2_lj.cpp
+++ b/src/d2_lj.cpp
@@ -19,6 +19,7 @@ void D2_LJ::calc_aed(
   aed_type &aed,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r2_inv = 1.0/rij_sq;
   double r6_inv = r2_inv*r2_inv*r2_inv;
   double fc_ij = fcut->calc(rij);
@@ -34,6 +35,7 @@ void D2_LJ::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
 
   double r2_inv = 1.0/rij_sq;
   double r6_inv = r2_inv*r2_inv*r2_inv;
@@ -52,6 +54,7 @@ void D2_LJ::calc_all(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r2_inv = 1.0/rij_sq;
   double r6_inv = r2_inv*r2_inv*r2_inv;
 
diff --git a/src/d2_mie.cpp b/src/d2_mie.cpp
index 318b4dc..2d35258 100644
--- a/src/d2_mie.cpp
+++ b/src/d2_mie.cpp
@@ -24,6 +24,7 @@ void D2_MIE::calc_aed(
   aed_type &aed,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double rn_inv = 1.0/pow(rij,n);
   double rm_inv = 1.0/pow(rij,m);
   double fc_ij = fcut->calc(rij);
@@ -39,6 +40,7 @@ void D2_MIE::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double fc_ij = fcut->calc(rij);
   double fcp_ij = fcut->calc_prime(rij);
 
@@ -54,6 +56,7 @@ void D2_MIE::calc_all(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
 
   double rn_inv = 1.0/pow(rij,n);
   double rm_inv = 1.0/pow(rij,m);
diff --git a/src/d2_mjoin.cpp b/src/d2_mjoin.cpp
index 44edea5..7988ea3 100644
--- a/src/d2_mjoin.cpp
+++ b/src/d2_mjoin.cpp
@@ -22,8 +22,6 @@ D2_mJoin::D2_mJoin(Config &c) : D2_Base(c) {
 void D2_mJoin::calc_aed(const int Zi, const int Zj, const double rij,
                         const double rij_sq, aed_type &aed, const double scale) {
   for (auto d : ds) {
-    if (!d->is_init_for_atoms(Zi,Zj)) continue;
-    if (rij > d->get_rcut()) continue;
     d->calc_aed(Zi, Zj, rij, rij_sq, aed, scale);
   }
 }
@@ -31,8 +29,6 @@ void D2_mJoin::calc_aed(const int Zi, const int Zj, const double rij,
 void D2_mJoin::calc_dXijdri(const int Zi, const int Zj, const double rij,
                             const double rij_sq, fd_type &fd_ij, const double scale) {
   for (auto d : ds) {
-    if (!d->is_init_for_atoms(Zi,Zj)) continue;
-    if (rij > d->get_rcut()) continue;
     d->calc_dXijdri(Zi, Zj, rij, rij_sq, fd_ij, scale);
   }
 }
@@ -40,8 +36,6 @@ void D2_mJoin::calc_dXijdri(const int Zi, const int Zj, const double rij,
 void D2_mJoin::calc_all(const int Zi, const int Zj, const double rij,
                         const double rij_sq, aed_type &aed, fd_type &fd_ij, const double scale) {
   for (auto d : ds) {
-    if (!d->is_init_for_atoms(Zi,Zj)) continue;
-    if (rij > d->get_rcut()) continue;
     d->calc_all(Zi, Zj, rij, rij_sq, aed, fd_ij, scale);
   }
 }
diff --git a/src/d2_zbl.cpp b/src/d2_zbl.cpp
index 5c2e7fb..e158b5e 100644
--- a/src/d2_zbl.cpp
+++ b/src/d2_zbl.cpp
@@ -56,6 +56,7 @@ void D2_ZBL::calc_aed(
   aed_type &aed,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double fc_ij = fcut->calc(rij);
   double a_ = a[types[Zi]][types[Zj]];
   double x = rij / a_;
@@ -70,6 +71,7 @@ void D2_ZBL::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r_inv = 1.0/rij;
   double a_ = a[types[Zi]][types[Zj]];
   double x = rij / a_;
@@ -91,6 +93,7 @@ void D2_ZBL::calc_all(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r_inv = 1.0/rij;
   double a_ = a[types[Zi]][types[Zj]];
   double x = rij / a_;
diff --git a/src/dm_blip.cpp b/src/dm_blip.cpp
index f75dc44..c68b8ec 100644
--- a/src/dm_blip.cpp
+++ b/src/dm_blip.cpp
@@ -71,6 +71,8 @@ void DM_Blip::calc_dXijdri_dXjidri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
+
   const double x = vec_ij[0];
   const double y = vec_ij[1];
   const double z = vec_ij[2];
@@ -163,6 +165,7 @@ void DM_Blip::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   const double x = vec_ij[0];
   const double y = vec_ij[1];
   const double z = vec_ij[2];
@@ -243,6 +246,7 @@ void DM_Blip::calc_rho(
   rho_type& rhoi,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t ii=rfidx;
   double fc_ij = scale*weights[Zj]*fcut->calc(rij);
   for (int L=0; L<=Lmax; ++L) {
diff --git a/src/dm_ead.cpp b/src/dm_ead.cpp
index c241da9..b7419b3 100644
--- a/src/dm_ead.cpp
+++ b/src/dm_ead.cpp
@@ -78,6 +78,7 @@ void DM_EAD::calc_dXijdri_dXjidri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   const double x = vec_ij[0];
   const double y = vec_ij[1];
   const double z = vec_ij[2];
@@ -171,6 +172,7 @@ void DM_EAD::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   const double x = vec_ij[0];
   const double y = vec_ij[1];
   const double z = vec_ij[2];
@@ -251,6 +253,7 @@ void DM_EAD::calc_rho(
   rho_type& rhoi,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   size_t ii=rfidx;
   double fc_ij = scale*weights[Zj]*fcut->calc(rij);
   for (int L=0; L<=Lmax; ++L) {
diff --git a/src/dm_eam.cpp b/src/dm_eam.cpp
index 94c0c48..8450987 100644
--- a/src/dm_eam.cpp
+++ b/src/dm_eam.cpp
@@ -58,8 +58,8 @@ void DM_EAM::calc_aed(
 
 }
 void DM_EAM::calc_dXijdri_dXjidri(
-  const int ,
-  const int ,
+  const int Zi,
+  const int Zj,
   const double rij,
   const double,
   const Vec3d &vec_ij,
@@ -68,6 +68,7 @@ void DM_EAM::calc_dXijdri_dXjidri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r = rij;
   double rinv = 1/r;
   double p = r*ef.rdr + 1.0;
@@ -88,8 +89,8 @@ void DM_EAM::calc_dXijdri_dXjidri(
   fd_ij(fidx,2) = f*vec_ij[2]; 
 }
 void DM_EAM::calc_dXijdri(
-  const int ,
-  const int ,
+  const int Zi,
+  const int Zj,
   const double rij,
   const double,
   const Vec3d &vec_ij,
@@ -97,6 +98,7 @@ void DM_EAM::calc_dXijdri(
   fd_type &fd_ij,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r = rij;
   double rinv = 1/r;
   double p = r*ef.rdr + 1.0;
@@ -192,14 +194,15 @@ void DM_EAM::init_rhoi(rho_type& rhoi)
   rhoi.set_zero();
 }
 void DM_EAM::calc_rho(
-  const int ,
-  const int ,
+  const int Zi,
+  const int Zj,
   const double rij,
   const double,
   const Vec3d &,
   rho_type& rho,
   const double scale)
 {
+  if (!is_init_for_atoms(Zi,Zj) || rij > get_rcut()) return;
   double r = rij;
   double p = r*ef.rdr + 1.0;
   int m = static_cast<int> (p);
diff --git a/src/dm_mjoin.cpp b/src/dm_mjoin.cpp
index 7e6a402..f1d9cb0 100644
--- a/src/dm_mjoin.cpp
+++ b/src/dm_mjoin.cpp
@@ -37,9 +37,7 @@ void DM_mJoin::calc_dXijdri_dXjidri(
   fd_type &fd_ij,
   const double scale) {
   for (auto d : ds) {
-    if (d->is_init_for_atoms(Zi,Zj) && rij < d->get_rcut()) {
-      d->calc_dXijdri_dXjidri(Zi,Zj,rij,rij_sq,vec_ij,rhoi,rhoj,fd_ij,scale); 
-    }
+    d->calc_dXijdri_dXjidri(Zi,Zj,rij,rij_sq,vec_ij,rhoi,rhoj,fd_ij,scale); 
   }
 }
 void DM_mJoin::calc_dXijdri(
@@ -52,9 +50,7 @@ void DM_mJoin::calc_dXijdri(
   fd_type &fd_ij, 
   const double scale) {
   for (auto d : ds) { 
-    if (d->is_init_for_atoms(Zi,Zj) && rij < d->get_rcut()) {
-      d->calc_dXijdri(Zi,Zj,rij,rij_sq,vec_ij,rhoi,fd_ij,scale);
-    }
+    d->calc_dXijdri(Zi,Zj,rij,rij_sq,vec_ij,rhoi,fd_ij,scale);
   }
 }
 void DM_mJoin::init_rhoi(rho_type& rhoi) {
@@ -70,9 +66,7 @@ void DM_mJoin::calc_rho(
   rho_type& rho,
   const double scale) {
   for (auto d : ds) {
-    if (d->is_init_for_atoms(Zi,Zj) && rij < d->get_rcut()) {
       d->calc_rho(Zi,Zj,rij,rij_sq,vec_ij,rho,scale);
-    }
   }
 }
 std::string DM_mJoin::label() {
-- 
GitLab