From c2bcf791967727449d9a93e465df161a03c8e87e Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Sat, 18 Jan 2025 11:45:27 -0500
Subject: [PATCH] avoid divisions

---
 src/ML-POD/eapod.cpp | 20 ++++++++++++--------
 1 file changed, 12 insertions(+), 8 deletions(-)

diff --git a/src/ML-POD/eapod.cpp b/src/ML-POD/eapod.cpp
index 11dace1f28..0525cfdb55 100644
--- a/src/ML-POD/eapod.cpp
+++ b/src/ML-POD/eapod.cpp
@@ -982,10 +982,14 @@ double EAPOD::peratom_environment_descriptors(double *cb, double *bd, double *tm
     D[j] = 1.0 / sum;
   }
 
-  double sum = 0;
+  double sum = 0.0;
   for (int j = 0; j < nClusters; j++) sum += D[j];
   double sumD = sum;
-  for (int j = 0; j < nClusters; j++) P[j] = D[j]/sum;
+  if (sum != 0.0) {
+    for (int j = 0; j < nClusters; j++) P[j] = D[j]/sum;
+  } else {
+    for (int j = 0; j < nClusters; j++) P[j] = 0.0;
+  }
 
   int nc = nCoeffPerElement*(ti[0]-1);
   double ei = coeff[0 + nc];
@@ -1008,13 +1012,13 @@ double EAPOD::peratom_environment_descriptors(double *cb, double *bd, double *tm
   }
 
   for (int m = 0; m<Mdesc; m++) {
-    double S1 = 1/sumD;
-    double S2 = sumD*sumD;
+    double S1 = 1.0/sumD;
+    double S2 = S1*S1;
     double sum = 0.0;
     for (int j=0; j<nClusters; j++) {
       double dP_dB = 0.0;
       for (int k = 0; k < nClusters; k++) {
-        double dP_dD = -D[j] / S2;
+        double dP_dD = -D[j] * S2;
         if (k==j) dP_dD += S1;
         double dD_dB = 0.0;
         double D2 = 2 * D[k] * D[k];
@@ -2844,11 +2848,11 @@ void EAPOD::peratomenvironment_descriptors(double *P, double *dP_dR, double *B,
   DGEMM(&chn, &chn, &nClusters, &Mdesc, &nComponents, &alpha, dD_dpca, &nClusters, ProjMat, &nComponents, &beta, dD_dB, &nClusters);
 
   // calculate dP_dD
-  double S1 = 1 / sumD;
-  double S2 = sumD * sumD;
+  double S1 = 1.0 / sumD;
+  double S2 = S1 * S1;
   for (int k = 0; k < nClusters; k++) {
     for (int j = 0; j < nClusters; j++) {
-      dP_dD[j + k * nClusters] = -D[j] / S2;
+      dP_dD[j + k * nClusters] = -D[j] * S2;
     }
   }
   for (int j = 0; j < nClusters; j++) {