implementation of openACC expol energies

Author: momokchung

include/ff/hippo/expol.h

@@ -5,7 +5,8 @@
 namespace tinker {
 void expolData(RcOp);
 
-void alterpol();
+void alterpol(real (*polscale)[3][3], real (*polinv)[3][3]);
+void dexpol();
 
 enum class ExpolScr
 {

include/ff/hippomod.h

@@ -41,10 +41,6 @@ TINKER_EXTERN int nmdwexclude;
 TINKER_EXTERN int (*mdwexclude)[2];
 TINKER_EXTERN real (*mdwexclude_scale)[3];
 
-TINKER_EXTERN int ndexclude;
-TINKER_EXTERN int (*dexclude)[2];
-TINKER_EXTERN real* dexclude_scale;
-
 TINKER_EXTERN int nwexclude;
 TINKER_EXTERN int (*wexclude)[2];
 TINKER_EXTERN real* wexclude_scale;

include/seq/damp_hippo.h

@@ -1,4 +1,6 @@
 #pragma once
+#include "ff/hippo/expol.h"
+#include "ff/hippomod.h"
 #include "math/sinhc.h"
 #include "seq/seq.h"
 
@@ -608,4 +610,27 @@ inline void damp_rep(real* restrict dmpik, real r, real rr1, real r2, real rr3,
    }
    // clang-format on
 }
+
+#pragma acc routine seq
+SEQ_CUDA
+inline void damp_expl(
+   ExpolScr scrtyp, real& restrict s2, real& restrict ds2, real r, real sizik, real alphai, real alphak)
+{
+   real alphaik, dmpik2, dampik, dampik2, expik, s;
+
+   if (scrtyp == ExpolScr::S2U) {
+      alphaik = REAL_SQRT(alphai * alphak);
+      real inv2 = 1. / 2, inv3 = 1. / 3;
+      real one = 1.;
+      dmpik2 = inv2 * alphaik;
+      dampik = dmpik2 * r;
+      dampik2 = dampik * dampik;
+      expik = REAL_EXP(-dampik);
+      s = (one + dampik + dampik2 * inv3) * expik;
+      s2 = s * s;
+      ds2 = s * (-alphaik * inv3) * (dampik + dampik2) * expik;
+   }
+   s2 = sizik * s2;
+   ds2 = sizik * ds2;
+}
 }

include/seq/pair_alterpol.h

@@ -0,0 +1,47 @@
+#pragma once
+#include "math/switch.h"
+#include "seq/damp_hippo.h"
+#include "ff/hippomod.h"
+
+namespace tinker {
+#pragma acc routine seq
+SEQ_CUDA
+inline void pair_alterpol(ExpolScr scrtyp, real r, real r2, real pscale, real cut, real off, real xr, real yr, real zr,
+   real springi, real sizi, real alphai, real springk, real sizk, real alphak,
+   real ks2i[3][3], real ks2k[3][3])
+{
+   real cut2 = cut * cut;
+   real sizik = sizi * sizk;
+   real s2;
+   real ds2;
+
+   damp_expl(scrtyp, s2, ds2, r, sizik, alphai, alphak);
+
+   if (r2 > cut2) {
+      real taper, dtaper;
+      switchTaper5<0>(r, cut, off, taper, dtaper);
+      s2 = s2 * taper;
+   }
+   real p33i, p33k;
+   p33i = springi * s2 * pscale;
+   p33k = springk * s2 * pscale;
+
+   real ai[3], ak[3];
+
+   ai[0] = xr / r;
+   ai[1] = yr / r;
+   ai[2] = zr / r;
+
+   ak[0] = -ai[0];
+   ak[1] = -ai[1];
+   ak[2] = -ai[2];
+   #pragma acc loop seq
+   for (int i{0}; i < 3; ++i) {
+      #pragma acc loop seq
+      for (int j{0}; j < 3; ++j) {
+         ks2i[j][i] = p33i * ai[i] * ai[j];
+         ks2k[j][i] = p33k * ak[i] * ak[j];
+      }
+   }
+}
+}

src/acc/hippo/alterpol.cpp

@@ -1,78 +1,156 @@
-//#include "ff/amoebamod.h"
+#include "ff/amoebamod.h"
 #include "ff/atom.h"
 #include "ff/hippomod.h"
 #include "ff/image.h"
 #include "ff/nblist.h"
 //#include "ff/pme.h"
 #include "ff/switch.h"
-#include "seq/pair_polar_chgpen.h"
-//#include "tool/gpucard.h"
-#include <array>
-#include <fstream>
+#include "seq/pair_alterpol.h"
+#include <tinker/routines.h>
+// #include "seq/pair_polar_chgpen.h"
+#include "tool/gpucard.h"
+// #include <array>
+// #include <fstream>
 
 namespace tinker {
-void alterpol()
+void alterpol(real (*polscale)[3][3], real (*polinv)[3][3])
 {
    real cut = switchCut(Switch::REPULS);
    real off = switchOff(Switch::REPULS);
 
    const real off2 = off * off;
-   const int maxnlst = mlist_unit->maxnlst;   // question
-   const auto* mlst = mlist_unit.deviceptr(); // question
-
-   size_t bufsize = bufferSize(); // question
-   PairPolarGrad pgrad;           // question
+   const int maxnlst = mlist_unit->maxnlst;
+   const auto* mlst = mlist_unit.deviceptr();
 
    // initialize polscale and polinv
+   #pragma acc parallel loop independent async\
+               deviceptr(polscale)
    for (int i = 0; i < n; ++i) {
+      polscale[i][0][0] = 1.f;
+      polscale[i][0][1] = 0.f;
+      polscale[i][0][2] = 0.f;
+      polscale[i][1][0] = 0.f;
       polscale[i][1][1] = 1.f;
       polscale[i][1][2] = 0.f;
-      polscale[i][1][3] = 0.f;
+      polscale[i][2][0] = 0.f;
       polscale[i][2][1] = 0.f;
       polscale[i][2][2] = 1.f;
-      polscale[i][2][3] = 0.f;
-      polscale[i][3][1] = 0.f;
-      polscale[i][3][2] = 0.f;
-      polscale[i][3][3] = 1.f;
-      polinv[i][1][1] = 1.f;
-      polinv[i][1][2] = 0.f;
-      polinv[i][1][3] = 0.f;
-      polinv[i][2][1] = 0.f;
-      polinv[i][2][2] = 1.f;
-      polinv[i][2][3] = 0.f;
-      polinv[i][3][1] = 0.f;
-      polinv[i][3][2] = 0.f;
-      polinv[i][3][3] = 1.f;
    }
 
    // find variable polarizability scale matrix at each site
+   MAYBE_UNUSED int GRID_DIM = gpuGridSize(BLOCK_DIM);
+   #pragma acc parallel async num_gangs(GRID_DIM) vector_length(BLOCK_DIM)\
+               deviceptr(x,y,z,kpep,prepep,dmppep,lpep,mlst,polscale)
+   #pragma acc loop gang independent
    for (int i = 0; i < n; ++i) {
       real springi = kpep[i];
       real sizi = prepep[i];
       real alphai = dmppep[i];
       int epli = lpep[i];
+      real xi = x[i];
+      real yi = y[i];
+      real zi = z[i];
 
-      int nmlsti = mlst->nlst[i]; // question
-      int base = i * maxnlst;     // question
+      int nmlsti = mlst->nlst[i];
+      int base = i * maxnlst;
+      #pragma acc loop vector independent
       for (int kk = 0; kk < nmlsti; ++kk) {
-         int offset = kk & (bufsize - 1); // question (copied from epolar.cpp)
-         int k = mlst->lst[base + kk];    // question
+         int k = mlst->lst[base + kk];
+         real xr = x[k] - xi;
+         real yr = y[k] - yi;
+         real zr = z[k] - zi;
+         real r2 = image2(xr, yr, zr);
          int eplk = lpep[k];
-         if (epli || eplk) {
-            real xr = x[k] - x[i];
-            real yr = y[k] - y[i];
-            real zr = z[k] - z[i];
+         bool incl = (epli || eplk);
+         if (r2 <= off2 and incl) {
+            real r = REAL_SQRT(r2);
+            real springk = kpep[k];
+            real sizk = prepep[k];
+            real alphak = dmppep[k];
+            real ks2i[3][3], ks2k[3][3];
+            pair_alterpol(scrtyp, r, r2, 1, cut, off, xr, yr, zr, springi, sizi, alphai, springk,
+               sizk, alphak, ks2i, ks2k);
+            #pragma acc loop seq
+            for (int l = 0; l < 3; ++l) {
+               #pragma acc loop seq
+               for (int m = 0; m < 3; ++m) {
+                  polscale[i][m][l] = polscale[i][m][l] + ks2i[m][l];
+                  polscale[k][m][l] = polscale[k][m][l] + ks2k[m][l];
+               }
+            }
+         }
+      }
+   }
+
+   #pragma acc parallel loop independent async\
+               deviceptr(x,y,z,kpep,prepep,dmppep,lpep,mlst,mdwexclude,mdwexclude_scale,polscale)
+   for (int ii = 0; ii < nmdwexclude; ++ii) {
+      int i = mdwexclude[ii][0];
+      int k = mdwexclude[ii][1];
+      real dscale = mdwexclude_scale[ii][1] - 1;
+
+      real springi = kpep[i];
+      real sizi = prepep[i];
+      real alphai = dmppep[i];
+      int epli = lpep[i];
 
-            zero(pgrad); // question
-            real r2 = image2(xr, yr, zr);
-            if (r2 <= off2) {
-               real r = REAL_SQRT(r2);
-               real sizk = prepep[k];
-               real alphak = dmppep[k];
-               // call dampexpl
+      real springk = kpep[k];
+      real sizk = prepep[k];
+      real alphak = dmppep[k];
+      int eplk = lpep[k];
+
+      real xr = x[k] - x[i];
+      real yr = y[k] - y[i];
+      real zr = z[k] - z[i];
+
+      real r2 = image2(xr, yr, zr);
+      bool incl1 = dscale != 0;
+      bool incl2 = (epli || eplk);
+      if (r2 <= off2 and incl1 and incl2) {
+         real r = REAL_SQRT(r2);
+         real ks2i[3][3], ks2k[3][3];
+         pair_alterpol(scrtyp, r, r2, 1, cut, off, xr, yr, zr, springi, sizi, alphai, springk, sizk,
+            alphak, ks2i, ks2k);
+         #pragma acc loop seq
+         for (int l = 0; l < 3; ++l) {
+            #pragma acc loop seq
+            for (int m = 0; m < 3; ++m) {
+               polscale[i][m][l] = polscale[i][m][l] + ks2i[m][l];
+               polscale[k][m][l] = polscale[k][m][l] + ks2k[m][l];
             }
          }
       }
    }
+   
+   // invert
+   #pragma acc parallel loop independent async\
+               deviceptr(polscale,polinv)
+   for (int i = 0; i < n; ++i) {
+      real tmp[3][3];
+      #pragma acc loop seq
+      for (int j = 0; j < 3; ++j) {
+         #pragma acc loop seq
+         for (int k = 0; k < 3; ++k) {
+            tmp[j][k] = polscale[i][j][k];
+         }
+      }
+      real det;
+      det = tmp[0][0] * (tmp[1][1] * tmp[2][2] - tmp[1][2] * tmp[2][1]) -
+         tmp[1][0] * (tmp[0][1] * tmp[2][2] - tmp[2][1] * tmp[0][2]) +
+         tmp[2][0] * (tmp[0][1] * tmp[1][2] - tmp[1][1] * tmp[0][2]);
+      polinv[i][0][0] = (tmp[1][1] * tmp[2][2] - tmp[1][2] * tmp[2][1]) / det;
+      polinv[i][1][0] = (tmp[2][0] * tmp[1][2] - tmp[1][0] * tmp[2][2]) / det;
+      polinv[i][2][0] = (tmp[1][0] * tmp[2][1] - tmp[2][0] * tmp[1][1]) / det;
+      polinv[i][0][1] = (tmp[2][1] * tmp[0][2] - tmp[0][1] * tmp[2][2]) / det;
+      polinv[i][1][1] = (tmp[0][0] * tmp[2][2] - tmp[2][0] * tmp[0][2]) / det;
+      polinv[i][2][1] = (tmp[0][1] * tmp[2][0] - tmp[0][0] * tmp[2][1]) / det;
+      polinv[i][0][2] = (tmp[0][1] * tmp[1][2] - tmp[0][2] * tmp[1][1]) / det;
+      polinv[i][1][2] = (tmp[0][2] * tmp[1][0] - tmp[0][0] * tmp[1][2]) / det;
+      polinv[i][2][2] = (tmp[0][0] * tmp[1][1] - tmp[0][1] * tmp[1][0]) / det;
+   }
+}
+
+void dexpol() {
+
+}
 }
-}