exchpol openacc gradient and virial implementation

Author: momokchung

include/seq/pair_alterpol.h

@@ -47,6 +47,8 @@ inline void pair_alterpol(ExpolScr scrtyp, real r, real r2, real pscale, real cu
    }
 }
 
+#pragma acc routine seq
+SEQ_CUDA
 inline void pair_dexpol(ExpolScr scrtyp, real r, real r2, real pscale, real cut, real off, real xr,
    real yr, real zr, real uix, real uiy, real uiz, real ukx, real uky, real ukz, real springi,
    real sizi, real alphai, real springk, real sizk, real alphak, const real f, real frc[3])
@@ -62,8 +64,8 @@ inline void pair_dexpol(ExpolScr scrtyp, real r, real r2, real pscale, real cut,
    if (r2 > cut2) {
       real taper, dtaper;
       switchTaper5<1>(r, cut, off, taper, dtaper);
-      s2 = s2 * taper;
       ds2 = ds2 * taper + s2 * dtaper;
+      s2 = s2 * taper;
    }
    real s2i = springi * s2 * pscale;
    real s2k = springk * s2 * pscale;

src/acc/hippo/alterpol.cpp

@@ -105,8 +105,8 @@ void alterpol(real (*polscale)[3][3], real (*polinv)[3][3])
       if (r2 <= off2 and incl1 and incl2) {
          real r = REAL_SQRT(r2);
          real ks2i[3][3], ks2k[3][3];
-         pair_alterpol(scrtyp, r, r2, dscale, cut, off, xr, yr, zr, springi, sizi, alphai, springk, sizk,
-            alphak, ks2i, ks2k);
+         pair_alterpol(scrtyp, r, r2, dscale, 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
@@ -139,8 +139,8 @@ void alterpol(real (*polscale)[3][3], real (*polinv)[3][3])
    }
 }
 
-void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* depy, grad_prec* depz,
-   VirialBuffer restrict vir_ep)
+void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* depy,
+   grad_prec* depz, VirialBuffer restrict vir_ep)
 {
    auto do_v = vers & calc::virial;
    real cut = switchCut(Switch::REPULS);
@@ -154,11 +154,15 @@ void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* d
 
    const real f = 0.5f * electric / dielec;
 
+   MAYBE_UNUSED int GRID_DIM = gpuGridSize(BLOCK_DIM);
+   #pragma acc parallel async num_gangs(GRID_DIM) vector_length(BLOCK_DIM)\
+               deviceptr(x,y,z,polarity,kpep,prepep,dmppep,lpep,uind,depx,depy,depz,vir_ep,mlst,polscale)
+   #pragma acc loop gang independent
    for (int i = 0; i < n; ++i) {
       real xi = x[i];
       real yi = y[i];
       real zi = z[i];
-      real springi = kpep[i]/polarity[i];
+      real springi = kpep[i] / polarity[i];
       real sizi = prepep[i];
       real alphai = dmppep[i];
       int epli = lpep[i];
@@ -170,6 +174,7 @@ void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* d
 
       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);
          int k = mlst->lst[base + kk];
@@ -181,7 +186,7 @@ void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* d
          bool incl = (epli || eplk);
          if (r2 <= off2 and incl) {
             real r = REAL_SQRT(r2);
-            real springk = kpep[k]/polarity[k];
+            real springk = kpep[k] / polarity[k];
             real sizk = prepep[k];
             real alphak = dmppep[k];
             real ukx = uind[k][0];
@@ -213,6 +218,9 @@ void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* d
       atomic_add(gzi, depz, i);
    }
 
+   #pragma acc parallel loop independent async\
+               deviceptr(x,y,z,polarity,kpep,prepep,dmppep,lpep,uind,depx,depy,depz,\
+                         vir_ep,mlst,mdwexclude,mdwexclude_scale,polscale)
    for (int ii = 0; ii < nmdwexclude; ++ii) {
       int offset = ii & (bufsize - 1);
 
@@ -245,16 +253,16 @@ void dexpol(const int vers, const real (*uind)[3], grad_prec* depx, grad_prec* d
       if (r2 <= off2 and incl1 and incl2) {
          real r = REAL_SQRT(r2);
          real frc[3];
-         pair_dexpol(scrtyp, r, r2, dscale, cut, off, xr, yr, zr, uix, uiy, uiz, ukx, uky, ukz, springi,
-            sizi, alphai, springk, sizk, alphak, f, frc);
+         pair_dexpol(scrtyp, r, r2, dscale, cut, off, xr, yr, zr, uix, uiy, uiz, ukx, uky, ukz,
+            springi, sizi, alphai, springk, sizk, alphak, f, frc);
 
          atomic_add(frc[0], depx, i);
          atomic_add(frc[1], depy, i);
          atomic_add(frc[2], depz, i);
          atomic_add(-frc[0], depx, k);
          atomic_add(-frc[1], depy, k);
          atomic_add(-frc[2], depz, k);
-         
+
          if (do_v) {
             real vxx = -xr * frc[0];
             real vxy = -0.5f * (yr * frc[0] + xr * frc[1]);