Add chkpole in repel multipole +elam in chgtrn fix

Author: roseane-reis

src/acc/amoeba/rotpole.cpp

@@ -18,6 +18,13 @@ void rotpole_acc()
       rotpoleAtomI(i, rpole, pole, zaxis, x, y, z);
 }
 
+void chkrepole_acc()
+{
+   #pragma acc parallel loop independent async deviceptr(x,y,z,zaxis,repole)
+   for (int i = 0; i < n; ++i)
+      chkpoleAtomI(i, repole, zaxis, x, y, z);
+}
+
 void rotrepole_acc()
 {
    #pragma acc parallel loop independent async deviceptr(x,y,z,zaxis,rrepole,repole)

src/acc/aplus/echgtrn.cpp

@@ -10,7 +10,7 @@
 #include <cassert>
 
 namespace tinker {
-#define DEVICE_PTRS x, y, z, dectx, decty, dectz, chgct, dmpct, nct, ect, vir_ect
+#define DEVICE_PTRS x, y, z, dectx, decty, dectz, chgct, dmpct, mut, elam, nct, ect, vir_ect
 template <class Ver>
 static void echgtrnAplus_acc1()
 {
@@ -42,6 +42,7 @@ static void echgtrnAplus_acc1()
       real zi = z[i];
       real chgi = chgct[i];
       real alphai = dmpct[i];
+      int imut = mut[i];
 
       MAYBE_UNUSED real gxi = 0, gyi = 0, gzi = 0;
 
@@ -57,12 +58,14 @@ static void echgtrnAplus_acc1()
          real zr = z[k] - zi;
          real chgk = chgct[k];
          real alphak = dmpct[k];
+         int kmut = mut[k];
 
          real r2 = image2(xr, yr, zr);
          if (r2 <= off2) {
             real r = REAL_SQRT(r2);
+            real elambda = (imut || kmut ? elam : 1);
             MAYBE_UNUSED real e, de;
-            pair_chgtrn_aplus<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, e, de);
+            pair_chgtrn_aplus<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, elambda, e, de);
 
             if CONSTEXPR (do_a)
                if (e != 0)
@@ -121,8 +124,10 @@ static void echgtrnAplus_acc1()
       real zi = z[i];
       real chgi = chgct[i];
       real alphai = dmpct[i];
+      int imut = mut[i];
       real chgk = chgct[k];
       real alphak = dmpct[k];
+      int kmut = mut[k];
 
       real xr = x[k] - xi;
       real yr = y[k] - yi;
@@ -131,12 +136,13 @@ static void echgtrnAplus_acc1()
       real r2 = image2(xr, yr, zr);
       if (r2 <= off2) {
          real r = REAL_SQRT(r2);
+         real elambda = (imut || kmut ? elam : 1);
          MAYBE_UNUSED real e, de;
-         pair_chgtrn_aplus<do_g>(r, cut, off, mscale, f, alphai, chgi, alphak, chgk, e, de);
+         pair_chgtrn_aplus<do_g>(r, cut, off, mscale, f, alphai, chgi, alphak, chgk, elambda, e, de);
 
          if CONSTEXPR (do_a) {
             real e1, de1;
-            pair_chgtrn_aplus<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, e1, de1);
+            pair_chgtrn_aplus<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, elambda, e1, de1);
 
             if (mscale == -1 and e1 != 0)
                atomic_add(-1, nct, offset);

src/acc/hippo/echgtrn.cpp

@@ -10,7 +10,7 @@
 #include <cassert>
 
 namespace tinker {
-#define DEVICE_PTRS x, y, z, dectx, decty, dectz, chgct, dmpct, nct, ect, vir_ect
+#define DEVICE_PTRS x, y, z, dectx, decty, dectz, chgct, dmpct, mut, elam, nct, ect, vir_ect
 template <class Ver>
 static void echgtrn_acc1()
 {
@@ -42,6 +42,7 @@ static void echgtrn_acc1()
       real zi = z[i];
       real chgi = chgct[i];
       real alphai = dmpct[i];
+      int imut = mut[i];
 
       MAYBE_UNUSED real gxi = 0, gyi = 0, gzi = 0;
 
@@ -57,12 +58,14 @@ static void echgtrn_acc1()
          real zr = z[k] - zi;
          real chgk = chgct[k];
          real alphak = dmpct[k];
+         int kmut = mut[k];
 
          real r2 = image2(xr, yr, zr);
          if (r2 <= off2) {
             real r = REAL_SQRT(r2);
+            real elambda = (imut || kmut ? elam : 1);
             MAYBE_UNUSED real e, de;
-            pair_chgtrn<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, e, de);
+            pair_chgtrn<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, elambda, e, de);
 
             if CONSTEXPR (do_a)
                if (e != 0)
@@ -121,8 +124,10 @@ static void echgtrn_acc1()
       real zi = z[i];
       real chgi = chgct[i];
       real alphai = dmpct[i];
+      int imut = mut[k];
       real chgk = chgct[k];
       real alphak = dmpct[k];
+      int kmut = mut[k];
 
       real xr = x[k] - xi;
       real yr = y[k] - yi;
@@ -131,12 +136,13 @@ static void echgtrn_acc1()
       real r2 = image2(xr, yr, zr);
       if (r2 <= off2) {
          real r = REAL_SQRT(r2);
+         real elambda = (imut || kmut ? elam : 1);
          MAYBE_UNUSED real e, de;
-         pair_chgtrn<do_g>(r, cut, off, mscale, f, alphai, chgi, alphak, chgk, e, de);
+         pair_chgtrn<do_g>(r, cut, off, mscale, f, alphai, chgi, alphak, chgk, elambda, e, de);
 
          if CONSTEXPR (do_a) {
             real e1, de1;
-            pair_chgtrn<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, e1, de1);
+            pair_chgtrn<do_g>(r, cut, off, 1, f, alphai, chgi, alphak, chgk, elambda, e1, de1);
 
             if (mscale == -1 and e1 != 0)
                atomic_add(-1, nct, offset);

src/cu/amoeba/rotpole.cu

@@ -15,6 +15,19 @@ void chkpole_cu()
 {
    launch_k1s(g::s0, n, chkpole_cu1, n, pole, zaxis, x, y, z);
 }
+
+__global__
+void chkrepole_cu1(int n, real (*restrict repole)[MPL_TOTAL], LocalFrame* zaxis, const real* restrict x,
+   const real* restrict y, const real* restrict z)
+{
+   for (int i = ITHREAD; i < n; i += STRIDE)
+      chkpoleAtomI(i, repole, zaxis, x, y, z);
+}
+
+void chkrepole_cu()
+{
+   launch_k1s(g::s0, n, chkrepole_cu1, n, repole, zaxis, x, y, z);
+}
 }
 
 namespace tinker {

src/hippo/erepel.cpp

@@ -157,10 +157,10 @@ void erepelData(RcOp op)
    }
 }
 
-TINKER_FVOID2(acc1, cu1, chkpole);
-static void chkpole()
+TINKER_FVOID2(acc1, cu1, chkrepole);
+static void chkrepole()
 {
-   TINKER_FCALL2(acc1, cu1, chkpole);
+   TINKER_FCALL2(acc1, cu1, chkrepole);
 }
 
 TINKER_FVOID2(acc1, cu1, rotrepole);
@@ -176,7 +176,7 @@ void repoleInit(int vers)
    if (vers & calc::virial)
       darray::zero(g::q0, bufferSize(), vir_trq);
 
-   chkpole();
+   chkrepole();
    rotrepole();
 
 }