Skip to content

Commit e5ac642

Browse files
zhi-wangzhi-wang
authored
Merge pull request #213 from roseane-reis/hippo-dev
Support for lambda scaling in HIPPO
2 parents 2334c76 + 5df4b15 commit e5ac642

29 files changed

+450
-200
lines changed

‎ext/ext/y3/echgtrn_cu1.yaml‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,8 @@ EXTRA_PARAMS: |
2828
, real *restrict chgct
2929
, real *restrict dmpct
3030
, real f
31+
, const int* restrict mut
32+
, real elam
3133
# IMPLICIT_PARAMS:
3234

3335
I_VARIABLES:
@@ -36,12 +38,14 @@ I_VARIABLES:
3638
- register real zi from:z
3739
- register real chgi from:chgct
3840
- register real alphai from:dmpct
41+
- register int imut from:mut
3942
K_VARIABLES:
4043
- register real xk from:x
4144
- register real yk from:y
4245
- register real zk from:z
4346
- register real chgk from:chgct
4447
- register real alphak from:dmpct
48+
- register int kmut from:mut
4549
I_FORCE:
4650
- register real gxi addto:gx
4751
- register real gyi addto:gy
@@ -58,11 +62,12 @@ SCALED_PAIRWISE_INTERACTION: |
5862
real r2 = image2(xr, yr, zr);
5963
if (r2 <= off * off and incl) {
6064
real r = REAL_SQRT(r2);
65+
real elambda = (imut || kmut ? elam : 1);
6166
e_prec e, de;
6267
if CONSTEXPR (CT == Chgtrn::SEPARATE)
63-
pair_chgtrn<do_g> (r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, e, de);
68+
pair_chgtrn<do_g> (r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, elambda, e, de);
6469
else if CONSTEXPR (CT == Chgtrn::COMBINED)
65-
pair_chgtrn_aplus<do_g>(r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, e, de);
70+
pair_chgtrn_aplus<do_g>(r, cut, off, scalea, f, @alphai@, @chgi@, @alphak@, @chgk@, elambda, e, de);
6671
if CONSTEXPR (do_a)
6772
if (e != 0 and scalea != 0)
6873
nctl += 1;
@@ -98,11 +103,12 @@ FULL_PAIRWISE_INTERACTION: |
98103
real r2 = image2(xr, yr, zr);
99104
if (r2 <= off * off and incl) {
100105
real r = REAL_SQRT(r2);
106+
real elambda = (imut || kmut ? elam : 1);
101107
e_prec e, de;
102108
if CONSTEXPR (CT == Chgtrn::SEPARATE)
103-
pair_chgtrn<do_g> (r, cut, off, 1, f, @alphai@, @chgi@, @alphak@, @chgk@, e, de);
109+
pair_chgtrn<do_g> (r, cut, off, 1, f, @alphai@, @chgi@, @alphak@, @chgk@, elambda, e, de);
104110
else if CONSTEXPR (CT == Chgtrn::COMBINED)
105-
pair_chgtrn_aplus<do_g>(r, cut, off, 1, f, @alphai@, @chgi@, @alphak@, @chgk@, e, de);
111+
pair_chgtrn_aplus<do_g>(r, cut, off, 1, f, @alphai@, @chgi@, @alphak@, @chgk@, elambda, e, de);
106112
if CONSTEXPR (do_a)
107113
if (e != 0)
108114
nctl += 1;

‎ext/ext/y3/edisp_cu1.yaml‎

Lines changed: 19 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -26,20 +26,23 @@ SCALE_1X_TYPE: real_const_array
2626

2727
EXTRA_PARAMS: |
2828
, const real* restrict csix, const real* restrict adisp, real aewald
29+
, const int* restrict mut, real vlam, Vdw vcouple
2930
# IMPLICIT_PARAMS:
3031

3132
I_VARIABLES:
32-
- register real xi from:x
33-
- register real yi from:y
34-
- register real zi from:z
35-
- register real ci from:csix
36-
- register real ai from:adisp
33+
- register real xi from:x
34+
- register real yi from:y
35+
- register real zi from:z
36+
- register real ci from:csix
37+
- register real ai from:adisp
38+
- register int imut from:mut
3739
K_VARIABLES:
38-
- register real xk from:x
39-
- register real yk from:y
40-
- register real zk from:z
41-
- register real ck from:csix
42-
- register real ak from:adisp
40+
- register real xk from:x
41+
- register real yk from:y
42+
- register real zk from:z
43+
- register real ck from:csix
44+
- register real ak from:adisp
45+
- register int kmut from:mut
4346
I_FORCE:
4447
- register real gxi addto:gx
4548
- register real gyi addto:gy
@@ -58,8 +61,9 @@ SCALED_PAIRWISE_INTERACTION: |
5861
real r = REAL_SQRT(r2);
5962
real rr1 = REAL_RECIP(r);
6063
real e, de;
61-
pair_disp<do_g, DTYP, 0>(r, r2, rr1,
62-
scalea, aewald, @ci@, @ai@, @ck@, @ak@, cut, off, e, de);
64+
real vlambda = pair_vlambda(vlam, vcouple, @imut@, @kmut@);
65+
pair_disp<do_g, DTYP, 0, 1>(r, r2, rr1,
66+
scalea, aewald, @ci@, @ai@, @ck@, @ak@, vlambda, cut, off, e, de);
6367
if CONSTEXPR (do_e) {
6468
edtl += floatTo<ebuf_prec>(e);
6569
if CONSTEXPR (do_a) {
@@ -98,8 +102,9 @@ FULL_PAIRWISE_INTERACTION: |
98102
real r = REAL_SQRT(r2);
99103
real rr1 = REAL_RECIP(r);
100104
real e, de;
101-
pair_disp<do_g, DTYP, 1>(r, r2, rr1,
102-
1, aewald, @ci@, @ai@, @ck@, @ak@, cut, off, e, de);
105+
real vlambda = pair_vlambda(vlam, vcouple, @imut@, @kmut@);
106+
pair_disp<do_g, DTYP, 1, 1>(r, r2, rr1,
107+
1, aewald, @ci@, @ai@, @ck@, @ak@, vlambda, cut, off, e, de);
103108
if CONSTEXPR (do_e) {
104109
edtl += floatTo<ebuf_prec>(e);
105110
if CONSTEXPR (do_a) {

‎ext/ext/y3/erepel_cu1.yaml‎

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,12 @@ EXTRA_PARAMS: |
2929
, real *restrict trqy
3030
, real *restrict trqz
3131
, const real(*restrict rpole)[10]
32-
, const real *restrict sizpr
33-
, const real *restrict elepr
32+
, const real* restrict sizpr
33+
, const real* restrict elepr
3434
, const real* restrict dmppr
35+
, const int* restrict mut
36+
, real vlam
37+
, Vdw vcouple
3538
# IMPLICIT_PARAMS:
3639

3740
I_VARIABLES:
@@ -51,6 +54,7 @@ I_VARIABLES:
5154
- shared real sizi from:sizpr
5255
- shared real dmpi from:dmppr
5356
- shared real vali from:elepr
57+
- shared int imut from:mut
5458
K_VARIABLES:
5559
- shared real xk from:x
5660
- shared real yk from:y
@@ -68,6 +72,7 @@ K_VARIABLES:
6872
- register real sizk from:sizpr
6973
- register real dmpk from:dmppr
7074
- register real valk from:elepr
75+
- register int kmut from:mut
7176
I_FORCE:
7277
- register real gxi addto:gx
7378
- register real gyi addto:gy
@@ -95,12 +100,13 @@ FULL_PAIRWISE_INTERACTION: |
95100
96101
real r2 = image2(xr, yr, zr);
97102
if (r2 <= off * off and incl) {
98-
pair_repel<do_g>( //
99-
r2, scalea, cut, off, xr, yr, zr,
103+
real vlambda = pair_vlambda(vlam, vcouple, @imut@, @kmut@);
104+
pair_repel<do_g, 1>( //
105+
r2, scalea, vlambda, cut, off, xr, yr, zr,
100106
@sizi@, @dmpi@, @vali@, @ci@, @dix@, @diy@, @diz@, @qixx@, @qixy@,
101107
@qixz@, @qiyy@, @qiyz@, @qizz@,
102108
@sizk@, @dmpk@, @valk@, @ck@, @dkx@, @dky@, @dkz@, @qkxx@, @qkxy@,
103-
@qkxz@, @qkyy@, @qkyz@, @qkzz@,
109+
@qkxz@, @qkyy@, @qkyz@, @qkzz@,
104110
e, pgrad);
105111
106112
if CONSTEXPR (do_a)

‎include/ff/elec.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -75,6 +75,7 @@ void elecData(RcOp);
7575

7676
TINKER_EXTERN real electric;
7777
TINKER_EXTERN real dielec;
78+
TINKER_EXTERN real elam;
7879

7980
/// \}
8081
}

‎include/ff/evdw.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,8 @@ enum class Vdw : int
5959

6060
namespace tinker {
6161
/// \ingroup vdw
62+
void vdwSoftcoreData(RcOp);
63+
/// \ingroup vdw
6264
void evdwData(RcOp);
6365
/// \ingroup vdw
6466
void evdw(int vers);

‎include/ff/hippo/echgtrn.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
#pragma once
2+
#include "ff/evdw.h"
23
#include "tool/rcman.h"
34

45
namespace tinker {

‎include/ff/hippo/edisp.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
#pragma once
22
#include "ff/energybuffer.h"
3+
#include "ff/evdw.h"
34
#include "tool/rcman.h"
45

56
namespace tinker {

‎include/ff/hippo/erepel.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,16 @@
11
#pragma once
2+
#include "ff/amoeba/mpole.h"
23
#include "ff/energybuffer.h"
4+
#include "ff/evdw.h"
35
#include "tool/rcman.h"
46

57
namespace tinker {
68
/// \ingroup repel
79
void erepelData(RcOp);
810
/// \ingroup repel
911
void erepel(int vers);
12+
/// \ingroup repel
13+
void repoleInit(int vers);
1014
}
1115

1216
//====================================================================//
@@ -16,6 +20,8 @@ void erepel(int vers);
1620
//====================================================================//
1721

1822
namespace tinker {
23+
TINKER_EXTERN real (*repole)[MPL_TOTAL];
24+
TINKER_EXTERN real (*rrepole)[MPL_TOTAL];
1925
TINKER_EXTERN real* sizpr;
2026
TINKER_EXTERN real* dmppr;
2127
TINKER_EXTERN real* elepr;

‎include/seq/pair_chgtrn.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,16 +8,16 @@ namespace tinker {
88
template <bool DO_G>
99
SEQ_CUDA
1010
void pair_chgtrn(real r, real cut, real off, real mscale, real f, real alphai,
11-
real chgi, real alphak, real chgk, e_prec& restrict e, e_prec& restrict de)
11+
real chgi, real alphak, real chgk, real elambda, e_prec& restrict e, e_prec& restrict de)
1212
{
1313
f *= mscale;
1414
real expi = REAL_EXP(-alphai * r);
1515
real expk = REAL_EXP(-alphak * r);
1616
e = -chgi * expk - chgk * expi;
17-
e *= f;
17+
e *= f * elambda;
1818
if CONSTEXPR (DO_G) {
1919
de = chgi * expk * alphak + chgk * expi * alphai;
20-
de *= f;
20+
de *= f * elambda;
2121
}
2222
if (r > cut) {
2323
real taper, dtaper;

‎include/seq/pair_disp.h‎

Lines changed: 60 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -4,19 +4,31 @@
44
#include "math/switch.h"
55
#include "seq/add.h"
66
#include "seq/damp_hippodisp.h"
7+
#include "seq/pair_vlambda.h"
78
#include "seq/seq.h"
89

910
namespace tinker {
1011
#pragma acc routine seq
1112
template <bool DO_G, class DTYP, int SCALE>
1213
SEQ_CUDA
13-
void pair_disp_obsolete(real r, real r2, real rr1, //
14-
real dspscale, real aewald, real ci, real ai, real ck, real ak, real edcut,
15-
real edoff, real& restrict e, real& restrict de)
14+
void pair_disp_obsolete(real r,
15+
real r2,
16+
real rr1,
17+
real dspscale,
18+
real aewald,
19+
real ci,
20+
real ai,
21+
real ck,
22+
real ak,
23+
real edcut,
24+
real edoff,
25+
real& restrict e,
26+
real& restrict de)
1627
{
1728
if (r > edoff) {
1829
e = 0;
19-
if CONSTEXPR (DO_G) de = 0;
30+
if CONSTEXPR (DO_G)
31+
de = 0;
2032
return;
2133
}
2234

@@ -62,10 +74,12 @@ void pair_disp_obsolete(real r, real r2, real rr1, //
6274
expi = REAL_EXP(-di);
6375
real term = ((((di + 5) * di + 17) * di / 96 + 0.5f) * di + 1) * di + 1;
6476
damp = 1 - term * expi;
65-
if CONSTEXPR (DO_G) ddamp = ai * expi * di2 * ((di2 - 3) * di - 3) / 96;
77+
if CONSTEXPR (DO_G)
78+
ddamp = ai * expi * di2 * ((di2 - 3) * di - 3) / 96;
6679
}
6780

68-
if CONSTEXPR (SCALE == 1) dspscale = 1;
81+
if CONSTEXPR (SCALE == 1)
82+
dspscale = 1;
6983

7084
if CONSTEXPR (eq<DTYP, DEWALD>()) {
7185
real ralpha2 = r2 * aewald * aewald;
@@ -75,8 +89,7 @@ void pair_disp_obsolete(real r, real r2, real rr1, //
7589
e = -ci * ck * rr6 * (dspscale * damp * damp + expa - 1);
7690
if CONSTEXPR (DO_G) {
7791
real rterm = -ralpha2 * ralpha2 * ralpha2 * rr1 * expterm;
78-
de = -6 * e * rr1
79-
- ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
92+
de = -6 * e * rr1 - ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
8093
}
8194
} else if CONSTEXPR (eq<DTYP, NON_EWALD_TAPER>()) {
8295
e = -ci * ck * rr6;
@@ -88,24 +101,38 @@ void pair_disp_obsolete(real r, real r2, real rr1, //
88101
if (r > edcut) {
89102
real taper, dtaper;
90103
switchTaper5<DO_G>(r, edcut, edoff, taper, dtaper);
91-
if CONSTEXPR (DO_G) de = e * dtaper + de * taper;
104+
if CONSTEXPR (DO_G)
105+
de = e * dtaper + de * taper;
92106
e = e * taper;
93107
}
94108
e *= dspscale;
95-
if CONSTEXPR (DO_G) de *= dspscale;
109+
if CONSTEXPR (DO_G)
110+
de *= dspscale;
96111
}
97112
}
98113

99114
#pragma acc routine seq
100-
template <bool DO_G, class DTYP, int SCALE>
115+
template <bool DO_G, class DTYP, int SCALE, int SOFTCORE>
101116
SEQ_CUDA
102-
void pair_disp(real r, real r2, real rr1, //
103-
real dspscale, real aewald, real ci, real ai, real ck, real ak, real edcut,
104-
real edoff, real& restrict e, real& restrict de)
117+
void pair_disp(real r,
118+
real r2,
119+
real rr1,
120+
real dspscale,
121+
real aewald,
122+
real ci,
123+
real ai,
124+
real ck,
125+
real ak,
126+
real vlambda,
127+
real edcut,
128+
real edoff,
129+
real& restrict e,
130+
real& restrict de)
105131
{
106132
if (r > edoff) {
107133
e = 0;
108-
if CONSTEXPR (DO_G) de = 0;
134+
if CONSTEXPR (DO_G)
135+
de = 0;
109136
return;
110137
}
111138

@@ -114,9 +141,20 @@ void pair_disp(real r, real r2, real rr1, //
114141
real dmpik[2], damp, ddamp;
115142
damp_hippodisp<DO_G>(dmpik, r, rr1, ai, ak);
116143
damp = dmpik[0];
117-
if CONSTEXPR (DO_G) ddamp = dmpik[1];
144+
if CONSTEXPR (DO_G)
145+
ddamp = dmpik[1];
146+
147+
if CONSTEXPR (SCALE == 1)
148+
dspscale = 1;
149+
150+
// set use of lambda scaling for decoupling or annihilation
151+
if CONSTEXPR (SOFTCORE) {
152+
real vlambda2 = vlambda * vlambda;
153+
real vlambda3 = vlambda2 * vlambda;
154+
real vterm = (vlambda2 * vlambda2) / REAL_SQRT(1.0 + vlambda2 - vlambda3);
155+
dspscale *= vterm;
156+
}
118157

119-
if CONSTEXPR (SCALE == 1) dspscale = 1;
120158
if CONSTEXPR (eq<DTYP, DEWALD>()) {
121159
real ralpha2 = r2 * aewald * aewald;
122160
real term = 1 + ralpha2 + 0.5f * ralpha2 * ralpha2;
@@ -125,8 +163,7 @@ void pair_disp(real r, real r2, real rr1, //
125163
e = -ci * ck * rr6 * (dspscale * damp * damp + expa - 1);
126164
if CONSTEXPR (DO_G) {
127165
real rterm = -ralpha2 * ralpha2 * ralpha2 * rr1 * expterm;
128-
de = -6 * e * rr1
129-
- ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
166+
de = -6 * e * rr1 - ci * ck * rr6 * (rterm + 2 * dspscale * damp * ddamp);
130167
}
131168
} else if CONSTEXPR (eq<DTYP, NON_EWALD_TAPER>()) {
132169
e = -ci * ck * rr6;
@@ -138,11 +175,13 @@ void pair_disp(real r, real r2, real rr1, //
138175
if (r > edcut) {
139176
real taper, dtaper;
140177
switchTaper5<DO_G>(r, edcut, edoff, taper, dtaper);
141-
if CONSTEXPR (DO_G) de = e * dtaper + de * taper;
178+
if CONSTEXPR (DO_G)
179+
de = e * dtaper + de * taper;
142180
e = e * taper;
143181
}
144182
e *= dspscale;
145-
if CONSTEXPR (DO_G) de *= dspscale;
183+
if CONSTEXPR (DO_G)
184+
de *= dspscale;
146185
}
147186
}
148187
}

0 commit comments

Comments
 (0)