Merge pull request #213 from roseane-reis/hippo-dev

Support for lambda scaling in HIPPO

ext/ext/y3/echgtrn_cu1.yaml

@@ -28,6 +28,8 @@ EXTRA_PARAMS: |
   , real *restrict chgct
   , real *restrict dmpct
   , real f
+  , const int* restrict mut
+  , real elam
 # IMPLICIT_PARAMS:
 
 I_VARIABLES:
@@ -36,12 +38,14 @@ I_VARIABLES:
   - register real zi     from:z
   - register real chgi   from:chgct
   - register real alphai from:dmpct
+  - register int  imut   from:mut
 K_VARIABLES:
   - register real xk     from:x
   - register real yk     from:y
   - register real zk     from:z
   - register real chgk   from:chgct
   - register real alphak from:dmpct
+  - register int  kmut   from:mut
 I_FORCE:
   - register real gxi addto:gx
   - register real gyi addto:gy
@@ -58,11 +62,12 @@ SCALED_PAIRWISE_INTERACTION: |
   real r2 = image2(xr, yr, zr);
   if (r2 <= off * off and incl) {
     real r = REAL_SQRT(r2);
+    real elambda = (imut || kmut ? elam : 1);
     e_prec e, de;
     if CONSTEXPR (CT == Chgtrn::SEPARATE)
-      pair_chgtrn<do_g>      (r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, e, de);
+      pair_chgtrn<do_g>      (r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, elambda, e, de);
     else if CONSTEXPR (CT == Chgtrn::COMBINED)
-      pair_chgtrn_aplus<do_g>(r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, e, de);
+      pair_chgtrn_aplus<do_g>(r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, elambda, e, de);
     if CONSTEXPR (do_a)
         if (e != 0 and scalea != 0)
           nctl += 1;
@@ -98,11 +103,12 @@ FULL_PAIRWISE_INTERACTION: |
   real r2 = image2(xr, yr, zr);
   if (r2 <= off * off and incl) {
     real r = REAL_SQRT(r2);
+    real elambda = (imut || kmut ? elam : 1);
     e_prec e, de;
     if CONSTEXPR (CT == Chgtrn::SEPARATE)
-      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);
     else if CONSTEXPR (CT == Chgtrn::COMBINED)
-      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)
           nctl += 1;

ext/ext/y3/edisp_cu1.yaml

@@ -26,20 +26,23 @@ SCALE_1X_TYPE: real_const_array
 
 EXTRA_PARAMS: |
   , const real* restrict csix, const real* restrict adisp, real aewald
+  , const int* restrict mut, real vlam, Vdw vcouple
 # IMPLICIT_PARAMS:
 
 I_VARIABLES:
-  - register real xi from:x
-  - register real yi from:y
-  - register real zi from:z
-  - register real ci from:csix
-  - register real ai from:adisp
+  - register real xi   from:x
+  - register real yi   from:y
+  - register real zi   from:z
+  - register real ci   from:csix
+  - register real ai   from:adisp
+  - register int  imut from:mut
 K_VARIABLES:
-  - register real xk from:x
-  - register real yk from:y
-  - register real zk from:z
-  - register real ck from:csix
-  - register real ak from:adisp
+  - register real xk   from:x
+  - register real yk   from:y
+  - register real zk   from:z
+  - register real ck   from:csix
+  - register real ak   from:adisp
+  - register int  kmut from:mut
 I_FORCE:
   - register real gxi addto:gx
   - register real gyi addto:gy
@@ -58,8 +61,9 @@ SCALED_PAIRWISE_INTERACTION: |
     real r = REAL_SQRT(r2);
     real rr1 = REAL_RECIP(r);
     real e, de;
-    pair_disp<do_g, DTYP, 0>(r, r2, rr1,
-      scalea, aewald, @ci@, @ai@, @ck@, @ak@, cut, off, e, de);
+    real vlambda = pair_vlambda(vlam, vcouple, @imut@, @kmut@);
+    pair_disp<do_g, DTYP, 0, 1>(r, r2, rr1,
+      scalea, aewald, @ci@, @ai@, @ck@, @ak@, vlambda, cut, off, e, de);
     if CONSTEXPR (do_e) {
       edtl += floatTo<ebuf_prec>(e);
       if CONSTEXPR (do_a) {
@@ -98,8 +102,9 @@ FULL_PAIRWISE_INTERACTION: |
     real r = REAL_SQRT(r2);
     real rr1 = REAL_RECIP(r);
     real e, de;
-    pair_disp<do_g, DTYP, 1>(r, r2, rr1,
-      1, aewald, @ci@, @ai@, @ck@, @ak@, cut, off, e, de);
+    real vlambda = pair_vlambda(vlam, vcouple, @imut@, @kmut@);
+    pair_disp<do_g, DTYP, 1, 1>(r, r2, rr1,
+      1, aewald, @ci@, @ai@, @ck@, @ak@, vlambda, cut, off, e, de);
     if CONSTEXPR (do_e) {
       edtl += floatTo<ebuf_prec>(e);
       if CONSTEXPR (do_a) {

ext/ext/y3/erepel_cu1.yaml

@@ -29,9 +29,12 @@ EXTRA_PARAMS: |
   , real *restrict trqy
   , real *restrict trqz
   , const real(*restrict rpole)[10]
-  , const real *restrict sizpr
-  , const real *restrict elepr
+  , const real* restrict sizpr
+  , const real* restrict elepr
   , const real* restrict dmppr
+  , const int*  restrict mut
+  , real vlam
+  , Vdw vcouple
 # IMPLICIT_PARAMS:
 
 I_VARIABLES:
@@ -51,6 +54,7 @@ I_VARIABLES:
   - shared real sizi from:sizpr
   - shared real dmpi from:dmppr
   - shared real vali from:elepr
+  - shared int  imut from:mut
 K_VARIABLES:
   - shared   real xk  from:x
   - shared   real yk  from:y
@@ -68,6 +72,7 @@ K_VARIABLES:
   - register real sizk from:sizpr
   - register real dmpk from:dmppr
   - register real valk from:elepr
+  - register int  kmut from:mut
 I_FORCE:
   - register real gxi addto:gx
   - register real gyi addto:gy
@@ -95,12 +100,13 @@ FULL_PAIRWISE_INTERACTION: |
 
   real r2 = image2(xr, yr, zr);
   if (r2 <= off * off and incl) {
-    pair_repel<do_g>( //
-      r2, scalea, cut, off, xr, yr, zr,
+    real vlambda = pair_vlambda(vlam, vcouple, @imut@, @kmut@);
+    pair_repel<do_g, 1>( //
+      r2, scalea, vlambda, cut, off, xr, yr, zr,
       @sizi@, @dmpi@, @vali@, @ci@, @dix@, @diy@, @diz@, @qixx@, @qixy@,
       @qixz@, @qiyy@, @qiyz@, @qizz@,
       @sizk@, @dmpk@, @valk@, @ck@, @dkx@, @dky@, @dkz@, @qkxx@, @qkxy@,
-      @qkxz@, @qkyy@, @qkyz@, @qkzz@,
+      @qkxz@, @qkyy@, @qkyz@, @qkzz@, 
       e, pgrad);
 
     if CONSTEXPR (do_a)

include/ff/elec.h

@@ -75,6 +75,7 @@ void elecData(RcOp);
 
 TINKER_EXTERN real electric;
 TINKER_EXTERN real dielec;
+TINKER_EXTERN real elam;
 
 /// \}
 }

include/ff/evdw.h

@@ -59,6 +59,8 @@ enum class Vdw : int
 
 namespace tinker {
 /// \ingroup vdw
+void vdwSoftcoreData(RcOp);
+/// \ingroup vdw
 void evdwData(RcOp);
 /// \ingroup vdw
 void evdw(int vers);

include/ff/hippo/echgtrn.h

@@ -1,4 +1,5 @@
 #pragma once
+#include "ff/evdw.h"
 #include "tool/rcman.h"
 
 namespace tinker {

include/ff/hippo/edisp.h

@@ -1,5 +1,6 @@
 #pragma once
 #include "ff/energybuffer.h"
+#include "ff/evdw.h"
 #include "tool/rcman.h"
 
 namespace tinker {

include/ff/hippo/erepel.h

@@ -1,12 +1,16 @@
 #pragma once
+#include "ff/amoeba/mpole.h"
 #include "ff/energybuffer.h"
+#include "ff/evdw.h"
 #include "tool/rcman.h"
 
 namespace tinker {
 /// \ingroup repel
 void erepelData(RcOp);
 /// \ingroup repel
 void erepel(int vers);
+/// \ingroup repel
+void repoleInit(int vers);
 }
 
 //====================================================================//
@@ -16,6 +20,8 @@ void erepel(int vers);
 //====================================================================//
 
 namespace tinker {
+TINKER_EXTERN real (*repole)[MPL_TOTAL];
+TINKER_EXTERN real (*rrepole)[MPL_TOTAL];
 TINKER_EXTERN real* sizpr;
 TINKER_EXTERN real* dmppr;
 TINKER_EXTERN real* elepr;

include/seq/pair_chgtrn.h

@@ -8,16 +8,16 @@ namespace tinker {
 template <bool DO_G>
 SEQ_CUDA
 void pair_chgtrn(real r, real cut, real off, real mscale, real f, real alphai,
-   real chgi, real alphak, real chgk, e_prec& restrict e, e_prec& restrict de)
+   real chgi, real alphak, real chgk, real elambda, e_prec& restrict e, e_prec& restrict de)
 {
    f *= mscale;
    real expi = REAL_EXP(-alphai * r);
    real expk = REAL_EXP(-alphak * r);
    e = -chgi * expk - chgk * expi;
-   e *= f;
+   e *= f * elambda;
    if CONSTEXPR (DO_G) {
       de = chgi * expk * alphak + chgk * expi * alphai;
-      de *= f;
+      de *= f * elambda;
    }
    if (r > cut) {
       real taper, dtaper;

include/seq/pair_disp.h

@@ -4,19 +4,31 @@
 #include "math/switch.h"
 #include "seq/add.h"
 #include "seq/damp_hippodisp.h"
+#include "seq/pair_vlambda.h"
 #include "seq/seq.h"
 
 namespace tinker {
 #pragma acc routine seq
 template <bool DO_G, class DTYP, int SCALE>
 SEQ_CUDA
-void pair_disp_obsolete(real r, real r2, real rr1, //
-   real dspscale, real aewald, real ci, real ai, real ck, real ak, real edcut,
-   real edoff, real& restrict e, real& restrict de)
+void pair_disp_obsolete(real r,
+                        real r2,
+                        real rr1,
+                        real dspscale,
+                        real aewald,
+                        real ci,
+                        real ai,
+                        real ck,
+                        real ak,
+                        real edcut,
+                        real edoff,
+                        real& restrict e,
+                        real& restrict de)
 {
    if (r > edoff) {
       e = 0;
-      if CONSTEXPR (DO_G) de = 0;
+      if CONSTEXPR (DO_G)
+         de = 0;
       return;
    }
 
@@ -62,10 +74,12 @@ void pair_disp_obsolete(real r, real r2, real rr1, //
       expi = REAL_EXP(-di);
       real term = ((((di + 5) * di + 17) * di / 96 + 0.5f) * di + 1) * di + 1;
       damp = 1 - term * expi;
-      if CONSTEXPR (DO_G) ddamp = ai * expi * di2 * ((di2 - 3) * di - 3) / 96;
+      if CONSTEXPR (DO_G)
+         ddamp = ai * expi * di2 * ((di2 - 3) * di - 3) / 96;
    }
 
-   if CONSTEXPR (SCALE == 1) dspscale = 1;
+   if CONSTEXPR (SCALE == 1)
+      dspscale = 1;
 
    if CONSTEXPR (eq<DTYP, DEWALD>()) {
       real ralpha2 = r2 * aewald * aewald;
@@ -75,8 +89,7 @@ void pair_disp_obsolete(real r, real r2, real rr1, //
       e = -ci * ck * rr6 * (dspscale * damp * damp + expa - 1);
       if CONSTEXPR (DO_G) {
          real rterm = -ralpha2 * ralpha2 * ralpha2 * rr1 * expterm;
-         de = -6 * e * rr1
-            - ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
+         de = -6 * e * rr1 - ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
       }
    } else if CONSTEXPR (eq<DTYP, NON_EWALD_TAPER>()) {
       e = -ci * ck * rr6;
@@ -88,24 +101,38 @@ void pair_disp_obsolete(real r, real r2, real rr1, //
       if (r > edcut) {
          real taper, dtaper;
          switchTaper5<DO_G>(r, edcut, edoff, taper, dtaper);
-         if CONSTEXPR (DO_G) de = e * dtaper + de * taper;
+         if CONSTEXPR (DO_G)
+            de = e * dtaper + de * taper;
          e = e * taper;
       }
       e *= dspscale;
-      if CONSTEXPR (DO_G) de *= dspscale;
+      if CONSTEXPR (DO_G)
+         de *= dspscale;
    }
 }
 
 #pragma acc routine seq
-template <bool DO_G, class DTYP, int SCALE>
+template <bool DO_G, class DTYP, int SCALE, int SOFTCORE>
 SEQ_CUDA
-void pair_disp(real r, real r2, real rr1, //
-   real dspscale, real aewald, real ci, real ai, real ck, real ak, real edcut,
-   real edoff, real& restrict e, real& restrict de)
+void pair_disp(real r,
+               real r2,
+               real rr1,
+               real dspscale,
+               real aewald,
+               real ci,
+               real ai,
+               real ck,
+               real ak,
+               real vlambda,
+               real edcut,
+               real edoff,
+               real& restrict e,
+               real& restrict de)
 {
    if (r > edoff) {
       e = 0;
-      if CONSTEXPR (DO_G) de = 0;
+      if CONSTEXPR (DO_G)
+         de = 0;
       return;
    }
 
@@ -114,9 +141,20 @@ void pair_disp(real r, real r2, real rr1, //
    real dmpik[2], damp, ddamp;
    damp_hippodisp<DO_G>(dmpik, r, rr1, ai, ak);
    damp = dmpik[0];
-   if CONSTEXPR (DO_G) ddamp = dmpik[1];
+   if CONSTEXPR (DO_G)
+      ddamp = dmpik[1];
+
+   if CONSTEXPR (SCALE == 1)
+      dspscale = 1;
+
+   // set use of lambda scaling for decoupling or annihilation
+   if CONSTEXPR (SOFTCORE) {
+      real vlambda2 = vlambda * vlambda;
+      real vlambda3 = vlambda2 * vlambda;
+      real vterm = (vlambda2 * vlambda2) / REAL_SQRT(1.0 + vlambda2 - vlambda3);
+      dspscale *= vterm;
+   }
 
-   if CONSTEXPR (SCALE == 1) dspscale = 1;
    if CONSTEXPR (eq<DTYP, DEWALD>()) {
       real ralpha2 = r2 * aewald * aewald;
       real term = 1 + ralpha2 + 0.5f * ralpha2 * ralpha2;
@@ -125,8 +163,7 @@ void pair_disp(real r, real r2, real rr1, //
       e = -ci * ck * rr6 * (dspscale * damp * damp + expa - 1);
       if CONSTEXPR (DO_G) {
          real rterm = -ralpha2 * ralpha2 * ralpha2 * rr1 * expterm;
-         de = -6 * e * rr1
-            - ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
+         de = -6 * e * rr1 - ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
       }
    } else if CONSTEXPR (eq<DTYP, NON_EWALD_TAPER>()) {
       e = -ci * ck * rr6;
@@ -138,11 +175,13 @@ void pair_disp(real r, real r2, real rr1, //
       if (r > edcut) {
          real taper, dtaper;
          switchTaper5<DO_G>(r, edcut, edoff, taper, dtaper);
-         if CONSTEXPR (DO_G) de = e * dtaper + de * taper;
+         if CONSTEXPR (DO_G)
+            de = e * dtaper + de * taper;
          e = e * taper;
       }
       e *= dspscale;
-      if CONSTEXPR (DO_G) de *= dspscale;
+      if CONSTEXPR (DO_G)
+         de *= dspscale;
    }
 }
 }