Merge pull request #203 from momokchung/exchpol

Exchange Polarization

.clang-format

@@ -10,6 +10,7 @@ BreakConstructorInitializers:          BeforeComma
 PointerAlignment:                      Left
 AllowShortBlocksOnASingleLine:         Always
 AllowShortFunctionsOnASingleLine:      Empty
+AllowShortIfStatementsOnASingleLine:   WithoutElse
 
 MaxEmptyLinesToKeep:                   1
 AccessModifierOffset:                  -3

ext/ext/yaml/alterpol_cu1.yaml

@@ -0,0 +1,89 @@
+KERNEL_NAME: alterpol_cu1
+
+CUT_DISTANCE:
+  - 'cut'
+
+OFF_DISTANCE:
+  - 'off'
+
+EXCLUDE_INFO:
+  - 'dinfo'
+
+SCALE_1X_TYPE: real3_const_array,1
+
+EXTRA_PARAMS: |
+  , real (*restrict polscale)[9]
+  , const real* restrict kpep, const real* restrict prepep
+  , const real* restrict dmppep, const int* restrict lpep
+  , ExpolScr scrtyp
+
+I_POSITION:
+  - def: shared real xi from:x
+  - def: shared real yi from:y
+  - def: shared real zi from:z
+K_POSITION:
+  - def: register real xk from:x
+  - def: register real yk from:y
+  - def: register real zk from:z
+
+I_FORCE:
+  - def: shared real psci00 addto:polscale,0
+  - def: shared real psci01 addto:polscale,1
+  - def: shared real psci02 addto:polscale,2
+  - def: shared real psci10 addto:polscale,3
+  - def: shared real psci11 addto:polscale,4
+  - def: shared real psci12 addto:polscale,5
+  - def: shared real psci20 addto:polscale,6
+  - def: shared real psci21 addto:polscale,7
+  - def: shared real psci22 addto:polscale,8
+K_FORCE:
+  - def: register real psck00 addto:polscale,0
+  - def: register real psck01 addto:polscale,1
+  - def: register real psck02 addto:polscale,2
+  - def: register real psck10 addto:polscale,3
+  - def: register real psck11 addto:polscale,4
+  - def: register real psck12 addto:polscale,5
+  - def: register real psck20 addto:polscale,6
+  - def: register real psck21 addto:polscale,7
+  - def: register real psck22 addto:polscale,8
+
+I_VARIABLES:
+  - def: shared real springi from:kpep
+  - def: shared real sizi from:prepep
+  - def: shared real alphai from:dmppep
+  - def: shared int epli from:lpep
+K_VARIABLES:
+  - def: register real springk from:kpep
+  - def: register real sizk from:prepep
+  - def: register real alphak from:dmppep
+  - def: register int eplk from:lpep
+
+FULL_PAIRWISE_INTERACTION: |
+  real xr = xk - @xi@;
+  real yr = yk - @yi@;
+  real zr = zk - @zi@;
+  real r2 = image2(xr, yr, zr);
+  if ((eplk or @epli@) and r2 <= off * off and incl) {
+    real r = REAL_SQRT(r2);
+    real ks2i[3][3], ks2k[3][3];
+    pair_alterpol(scrtyp, r, scaleb, cut, off, xr, yr, zr, @springi@, @sizi@, @alphai@,
+      springk, sizk, alphak, ks2i, ks2k);
+    @psci00@ += ks2i[0][0];
+    @psci01@ += ks2i[0][1];
+    @psci02@ += ks2i[0][2];
+    @psci10@ += ks2i[1][0];
+    @psci11@ += ks2i[1][1];
+    @psci12@ += ks2i[1][2];
+    @psci20@ += ks2i[2][0];
+    @psci21@ += ks2i[2][1];
+    @psci22@ += ks2i[2][2];
+    psck00 += ks2k[0][0];
+    psck01 += ks2k[0][1];
+    psck02 += ks2k[0][2];
+    psck10 += ks2k[1][0];
+    psck11 += ks2k[1][1];
+    psck12 += ks2k[1][2];
+    psck20 += ks2k[2][0];
+    psck21 += ks2k[2][1];
+    psck22 += ks2k[2][2];
+  }

ext/ext/yaml/dexpol.yaml

@@ -0,0 +1,100 @@
+KERNEL_NAME: dexpol_cu1
+TEMPLATE_PARAMS: template<class Ver>
+CONSTEXPR_FLAGS: |
+  constexpr bool do_v = Ver::v;
+
+VIRIAL:
+  - vep
+GRADIENT:
+  - gx
+  - gy
+  - gz
+
+CUT_DISTANCE:
+  - 'cut'
+
+OFF_DISTANCE:
+  - 'off'
+
+EXCLUDE_INFO:
+  - 'dinfo'
+
+SCALE_1X_TYPE: real3_const_array,1
+
+EXTRA_PARAMS: |
+  , const real* restrict polarity, const real (*restrict uind)[3]
+  , const real* restrict kpep, const real* restrict prepep
+  , const real* restrict dmppep, const int* restrict lpep
+  , ExpolScr scrtyp, real f
+
+I_POSITION:
+  - def: shared real xi from:x
+  - def: shared real yi from:y
+  - def: shared real zi from:z
+K_POSITION:
+  - def: register real xk from:x
+  - def: register real yk from:y
+  - def: register real zk from:z
+
+I_FORCE:
+  - def: shared real frcxi addto:gx
+  - def: shared real frcyi addto:gy
+  - def: shared real frczi addto:gz
+
+K_FORCE:
+  - def: register real frcxk addto:gx
+  - def: register real frcyk addto:gy
+  - def: register real frczk addto:gz
+
+I_VARIABLES:
+  - def: shared real uix from:uind,0
+  - def: shared real uiy from:uind,1
+  - def: shared real uiz from:uind,2
+  - def: shared real springi from:kpep
+  - def: shared real sizi from:prepep
+  - def: shared real alphai from:dmppep
+  - def: shared int epli from:lpep
+  - def: shared real poli from:polarity
+
+K_VARIABLES:
+  - def: register real ukx from:uind,0
+  - def: register real uky from:uind,1
+  - def: register real ukz from:uind,2
+  - def: register real springk from:kpep
+  - def: register real sizk from:prepep
+  - def: register real alphak from:dmppep
+  - def: register int eplk from:lpep
+  - def: register real polk from:polarity
+
+FULL_PAIRWISE_INTERACTION: |
+  real xr = xk - @xi@;
+  real yr = yk - @yi@;
+  real zr = zk - @zi@;
+  real r2 = image2(xr, yr, zr);
+  if ((eplk or @epli@) and r2 <= off * off and incl) {
+    real r = REAL_SQRT(r2);
+    real frc[3];
+    pair_dexpol(scrtyp, r, scaleb, cut, off, xr, yr, zr, @uix@, @uiy@, @uiz@, ukx, uky, ukz,
+               @springi@/@poli@, @sizi@, @alphai@, springk/polk, sizk, alphak, f, frc);
+    @frcxi@ += frc[0];
+    @frcyi@ += frc[1];
+    @frczi@ += frc[2];
+    frcxk -= frc[0];
+    frcyk -= frc[1];
+    frczk -= frc[2];
+
+    if CONSTEXPR (do_v) {
+      real vxx = -xr * frc[0];
+      real vxy = -0.5f * (yr * frc[0] + xr * frc[1]);
+      real vxz = -0.5f * (zr * frc[0] + xr * frc[2]);
+      real vyy = -yr * frc[1];
+      real vyz = -0.5f * (zr * frc[1] + yr * frc[2]);
+      real vzz = -zr * frc[2];
+      veptlxx += floatTo<vbuf_prec>(vxx);
+      veptlyx += floatTo<vbuf_prec>(vxy);
+      veptlzx += floatTo<vbuf_prec>(vxz);
+      veptlyy += floatTo<vbuf_prec>(vyy);
+      veptlzy += floatTo<vbuf_prec>(vyz);
+      veptlzz += floatTo<vbuf_prec>(vzz);
+    }
+  }

include/ff/cuinduce.h

@@ -0,0 +1,30 @@
+#pragma once
+#include "ff/precision.h"
+
+namespace tinker {
+// udir = polarity * field
+
+__global__
+void pcgUdirV1(int n, const real* polarity, //
+   real (*udir)[3], const real (*field)[3]);
+
+__global__
+void pcgUdirV2(int n, const real* polarity, //
+   real (*udir)[3], real (*udirp)[3], const real (*field)[3], const real (*fieldp)[3]);
+
+// r(0) = E - (1/polarity + Tu) u(0) = (udir - u(0))/polarity + mutual field
+
+__global__
+void pcgRsd0V1(int n, const real* polarity_inv, real (*rsd)[3], //
+   const real (*udir)[3], const real (*uind)[3], const real (*field)[3]);
+
+__global__
+void pcgRsd0V2(int n, const real* polarity_inv, real (*rsd)[3], real (*rsp)[3], //
+   const real (*udir)[3], const real (*udip)[3], const real (*uind)[3], const real (*uinp)[3],
+   const real (*field)[3], const real (*fielp)[3]);
+
+__global__
+void pcgRsd0V3(int n, const real* polarity_inv, real (*rsd)[3], //
+   const real (*udir)[3], const real (*uind)[3], const real (*field)[3],
+   const real (*polscale)[3][3]);
+}

include/ff/hippo/expol.h

@@ -0,0 +1,10 @@
+#pragma once
+#include "ff/energybuffer.h"
+#include "tool/rcman.h"
+
+namespace tinker {
+void expolData(RcOp);
+
+void alterpol(real (*polscale)[3][3], real (*polinv)[3][3]);
+void dexpol(int vers);
+}

include/ff/hippo/expolscr.h

@@ -0,0 +1,11 @@
+#pragma once
+
+namespace tinker {
+enum class ExpolScr
+{
+   NONE,
+   S2U,
+   S2,
+   G,
+};
+}

include/ff/hippomod.h

@@ -2,6 +2,7 @@
 #include "ff/energybuffer.h"
 #include "ff/hippo/chgpen.h"
 #include "ff/hippo/echgtrn.h"
+#include "ff/hippo/expolscr.h"
 
 // mplpot
 namespace tinker {
@@ -40,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;
@@ -65,3 +62,14 @@ TINKER_EXTERN virial_prec virial_ect[9];
 
 TINKER_EXTERN Chgtrn ctrntyp;
 }
+
+// expol
+namespace tinker {
+TINKER_EXTERN real* kpep;
+TINKER_EXTERN real* prepep;
+TINKER_EXTERN real* dmppep;
+TINKER_EXTERN int* lpep;
+TINKER_EXTERN real (*polscale)[3][3];
+TINKER_EXTERN real (*polinv)[3][3];
+TINKER_EXTERN ExpolScr scrtyp;
+}