Fix ghost atom handling

Author: marvinfriede

include/dftd_ncoord.h

@@ -200,4 +200,4 @@ extern inline double log_cn_cut(const double cn_max, const double cn);
 
 extern inline double dlog_cn_cut(const double cn_max, const double cn);
 
-}; // namespace dftd4
+} // namespace dftd4

include/dftd_parameters.h

@@ -63,10 +63,10 @@ static const double zeff[119]{
  * These values are actually never used in the code. Only r4r2 is used.
  */
 static const double r2r4[119]{
-  0.0,             //  dummy
-  8.0589,  3.4698, //  H,He
+  0.0,                                                           //  dummy
+  8.0589,  3.4698,                                               //  H,He
   29.0974, 14.8517, 11.8799, 7.8715,  5.5588,  4.7566,  3.8025,
-  3.1036, //  Li-Ne
+  3.1036,                                                        //  Li-Ne
   26.1552, 17.2304, 17.7210, 12.7442, 9.5361,  8.1652,  6.7463,
   5.6004,                                                        //  Na-Ar
   29.2012, 22.3934,                                              //  K,Ca

src/damping/dftd_atm.cpp

@@ -301,15 +301,15 @@ int get_atm_dispersion_derivs(
         dgyjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * yjk;
         dgzjk = c9ijk * (-dang * fdmp + ang * dfdmp) / r2jk * zjk;
 
-        gradient(3 * iat + 0) += -dgxij - dgxik;
-        gradient(3 * iat + 1) += -dgyij - dgyik;
-        gradient(3 * iat + 2) += -dgzij - dgzik;
-        gradient(3 * jat + 0) += dgxij - dgxjk;
-        gradient(3 * jat + 1) += dgyij - dgyjk;
-        gradient(3 * jat + 2) += dgzij - dgzjk;
-        gradient(3 * kat + 0) += dgxik + dgxjk;
-        gradient(3 * kat + 1) += dgyik + dgyjk;
-        gradient(3 * kat + 2) += dgzik + dgzjk;
+        gradient(3 * ii + 0) += -dgxij - dgxik;
+        gradient(3 * ii + 1) += -dgyij - dgyik;
+        gradient(3 * ii + 2) += -dgzij - dgzik;
+        gradient(3 * jj + 0) += dgxij - dgxjk;
+        gradient(3 * jj + 1) += dgyij - dgyjk;
+        gradient(3 * jj + 2) += dgzij - dgzjk;
+        gradient(3 * kk + 0) += dgxik + dgxjk;
+        gradient(3 * kk + 1) += dgyik + dgyjk;
+        gradient(3 * kk + 2) += dgzik + dgzjk;
 
         dEdcn(ii) -= e * 0.5 * (dc6dcn(ii, jj) / c6ij + dc6dcn(ii, kk) / c6ik);
         dEdcn(jj) -= e * 0.5 * (dc6dcn(jj, ii) / c6ij + dc6dcn(jj, kk) / c6jk);

src/damping/dftd_rational.cpp

@@ -18,16 +18,15 @@
 
 #include <cmath>
 
+#include "damping/dftd_atm.h"
+#include "damping/dftd_rational.h"
 #include "dftd_cblas.h"
 #include "dftd_dispersion.h"
 #include "dftd_eeq.h"
 #include "dftd_geometry.h"
 #include "dftd_matrix.h"
 #include "dftd_ncoord.h"
 #include "dftd_parameters.h"
-#include "damping/dftd_atm.h"
-#include "damping/dftd_rational.h"
-
 
 namespace dftd4 {
 
@@ -119,8 +118,8 @@ int get_dispersion2_energy(
       }
 
       e = -c6ij * edisp * 0.5;
-      energy(iat) += e;
-      energy(jat) += e;
+      energy(ii) += e;
+      energy(jj) += e;
     }
   }
 
@@ -164,7 +163,7 @@ int get_dispersion2_derivs(
 
       r4r2ij = 3.0 * r4r2[izp] * r4r2[jzp];
       r0ij = par.a1 * sqrt(r4r2ij) + par.a2;
-      c6ij = c6(iat, jat);
+      c6ij = c6(ii, jj);
 
       t6 = fdmpr_bj(6, r, r0ij);
       t8 = fdmpr_bj(8, r, r0ij);
@@ -199,12 +198,12 @@ int get_dispersion2_derivs(
       dEdq(ii) -= dc6dq(ii, jj) * edisp;
       dEdq(jj) -= dc6dq(jj, ii) * edisp;
 
-      gradient(3 * iat) += dgx;
-      gradient(3 * iat + 1) += dgy;
-      gradient(3 * iat + 2) += dgz;
-      gradient(3 * jat) -= dgx;
-      gradient(3 * jat + 1) -= dgy;
-      gradient(3 * jat + 2) -= dgz;
+      gradient(3 * ii) += dgx;
+      gradient(3 * ii + 1) += dgy;
+      gradient(3 * ii + 2) += dgz;
+      gradient(3 * jj) -= dgx;
+      gradient(3 * jj + 1) -= dgy;
+      gradient(3 * jj + 2) -= dgz;
     }
   }
   return EXIT_SUCCESS;

src/dftd_dispersion.cpp

@@ -65,9 +65,9 @@ int get_dispersion(
   cn.NewVector(nat);
   q.NewVector(nat);
   if (lgrad) {
-    dcndr.NewMatrix(3 * mol.NAtoms, nat);
-    dqdr.NewMatrix(3 * mol.NAtoms, nat);
-    gradient.NewVector(3 * mol.NAtoms);
+    dcndr.NewMatrix(3 * nat, nat);
+    dqdr.NewMatrix(3 * nat, nat);
+    gradient.NewVector(3 * nat);
   }
 
   // calculate partial charges from EEQ model
@@ -138,8 +138,7 @@ int get_dispersion(
   );
   if (info != EXIT_SUCCESS) return info;
 
-  if (lgrad) { BLAS_Add_Mat_x_Vec(gradient, dqdr, dEdq, false, 1.0); }
-
+  if (lgrad) BLAS_Add_Mat_x_Vec(gradient, dqdr, dEdq, false, 1.0);
   dqdr.DelMat();
 
   // ----------------------------
@@ -224,9 +223,15 @@ int get_dispersion(
 
   // write to input gradient
   if (lgrad) {
-    for (int i = 0; i != 3 * mol.NAtoms; i++) {
-      GRAD[i] = gradient(i);
+    for (int i = 0, ii = 0; i != mol.NAtoms; i++) {
+      ii = realIdx(i);
+      if (ii < 0) continue;
+
+      GRAD[3 * i] = gradient(3 * ii);
+      GRAD[3 * i + 1] = gradient(3 * ii + 1);
+      GRAD[3 * i + 2] = gradient(3 * ii + 2);
     }
+
     gradient.DelVec();
   }