Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pair_field.h
1#pragma once
2#include "ff/elec.h"
3#include "math/realn.h"
4#include "seq/add.h"
5#include "seq/damp.h"
6
7namespace tinker {
8#pragma acc routine seq
9template <class ETYP>
11void pair_dfield(real r2, real xr, real yr, real zr, real dscale,
12 real pscale, //
13 real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz,
14 real qiyy, real qiyz, real qizz, real pdi,
15 real pti, //
16 real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz,
17 real qkyy, real qkyz, real qkzz, real pdk,
18 real ptk, //
19 real aewald, real3& restrict fid, real3& restrict fip, real3& restrict fkd,
20 real3& restrict fkp)
21{
22 real r = REAL_SQRT(r2);
23 real invr1 = REAL_RECIP(r);
24 real rr2 = invr1 * invr1;
25
26 real scale3, scale5, scale7;
27 damp_thole3(r, pdi, pti, pdk, ptk, scale3, scale5, scale7);
28
29 real bn[4];
30 if CONSTEXPR (eq<ETYP, EWALD>()) damp_ewald<4>(bn, r, invr1, rr2, aewald);
31 real rr1 = invr1;
32 real rr3 = rr1 * rr2;
33 real rr5 = 3 * rr1 * rr2 * rr2;
34 real rr7 = 15 * rr1 * rr2 * rr2 * rr2;
35
36 real dir = dix * xr + diy * yr + diz * zr;
37 real qix = qixx * xr + qixy * yr + qixz * zr;
38 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
39 real qiz = qixz * xr + qiyz * yr + qizz * zr;
40 real qir = qix * xr + qiy * yr + qiz * zr;
41 real dkr = dkx * xr + dky * yr + dkz * zr;
42 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
43 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
44 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
45 real qkr = qkx * xr + qky * yr + qkz * zr;
46
47 real3 dixyz = make_real3(dix, diy, diz);
48 real3 dkxyz = make_real3(dkx, dky, dkz);
49 real3 qixyz = make_real3(qix, qiy, qiz);
50 real3 qkxyz = make_real3(qkx, qky, qkz);
51 real3 dr = make_real3(xr, yr, zr);
52 real c1;
53 real3 inci, inck;
54
55 // d-field
56
57 if CONSTEXPR (eq<ETYP, EWALD>()) {
58 bn[1] -= (1 - scale3) * rr3;
59 bn[2] -= (1 - scale5) * rr5;
60 bn[3] -= (1 - scale7) * rr7;
61 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
62 bn[1] = dscale * scale3 * rr3;
63 bn[2] = dscale * scale5 * rr5;
64 bn[3] = dscale * scale7 * rr7;
65 }
66
67 c1 = -(bn[1] * ck - bn[2] * dkr + bn[3] * qkr);
68 inci = c1 * dr - bn[1] * dkxyz + 2 * bn[2] * qkxyz;
69 fid += inci;
70
71 c1 = (bn[1] * ci + bn[2] * dir + bn[3] * qir);
72 inck = c1 * dr - bn[1] * dixyz - 2 * bn[2] * qixyz;
73 fkd += inck;
74
75 // p-field
76
77 if CONSTEXPR (eq<ETYP, EWALD>()) {
78 fip += inci;
79 fkp += inck;
80 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
81 if (pscale == dscale) {
82 fip += inci;
83 fkp += inck;
84 } else {
85 bn[1] = pscale * scale3 * rr3;
86 bn[2] = pscale * scale5 * rr5;
87 bn[3] = pscale * scale7 * rr7;
88
89 c1 = -(bn[1] * ck - bn[2] * dkr + bn[3] * qkr);
90 fip += c1 * dr - bn[1] * dkxyz + 2 * bn[2] * qkxyz;
91
92 c1 = (bn[1] * ci + bn[2] * dir + bn[3] * qir);
93 fkp += c1 * dr - bn[1] * dixyz - 2 * bn[2] * qixyz;
94 }
95 }
96}
97
98#pragma acc routine seq
99template <class ETYP>
101void pair_dfield_v2(real r2, real xr, real yr, real zr, real dscale,
102 real pscale, real aewald, //
103 real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz,
104 real qiyy, real qiyz, real qizz, real pdi,
105 real pti, //
106 real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz,
107 real qkyy, real qkyz, real qkzz, real pdk,
108 real ptk, //
109 real& restrict fidx, real& restrict fidy, real& restrict fidz,
110 real& restrict fipx, real& restrict fipy, real& restrict fipz,
111 real& restrict fkdx, real& restrict fkdy, real& restrict fkdz,
112 real& restrict fkpx, real& restrict fkpy, real& restrict fkpz)
113{
114 real r = REAL_SQRT(r2);
115 real invr1 = REAL_RECIP(r);
116 real rr2 = invr1 * invr1;
117
118 real scale3, scale5, scale7;
119 damp_thole3(r, pdi, pti, pdk, ptk, scale3, scale5, scale7);
120
121 real bn[4];
122 if CONSTEXPR (eq<ETYP, EWALD>()) damp_ewald<4>(bn, r, invr1, rr2, aewald);
123 real rr1 = invr1;
124 real rr3 = rr1 * rr2;
125 real rr5 = 3 * rr1 * rr2 * rr2;
126 real rr7 = 15 * rr1 * rr2 * rr2 * rr2;
127
128 real dir = dix * xr + diy * yr + diz * zr;
129 real qix = qixx * xr + qixy * yr + qixz * zr;
130 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
131 real qiz = qixz * xr + qiyz * yr + qizz * zr;
132 real qir = qix * xr + qiy * yr + qiz * zr;
133 real dkr = dkx * xr + dky * yr + dkz * zr;
134 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
135 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
136 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
137 real qkr = qkx * xr + qky * yr + qkz * zr;
138
139 real3 dixyz = make_real3(dix, diy, diz);
140 real3 dkxyz = make_real3(dkx, dky, dkz);
141 real3 qixyz = make_real3(qix, qiy, qiz);
142 real3 qkxyz = make_real3(qkx, qky, qkz);
143 real3 dr = make_real3(xr, yr, zr);
144 real c1;
145 real3 inci, inck;
146
147 // d-field
148 real bn1, bn2, bn3;
149 if CONSTEXPR (eq<ETYP, EWALD>()) {
150 bn1 = bn[1];
151 bn2 = bn[2];
152 bn3 = bn[3];
153 bn[1] = bn1 - (1 - dscale * scale3) * rr3;
154 bn[2] = bn2 - (1 - dscale * scale5) * rr5;
155 bn[3] = bn3 - (1 - dscale * scale7) * rr7;
156 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
157 bn[1] = dscale * scale3 * rr3;
158 bn[2] = dscale * scale5 * rr5;
159 bn[3] = dscale * scale7 * rr7;
160 }
161
162 c1 = -(bn[1] * ck - bn[2] * dkr + bn[3] * qkr);
163 inci = c1 * dr - bn[1] * dkxyz + 2 * bn[2] * qkxyz;
164 fidx += inci.x;
165 fidy += inci.y;
166 fidz += inci.z;
167
168 c1 = (bn[1] * ci + bn[2] * dir + bn[3] * qir);
169 inck = c1 * dr - bn[1] * dixyz - 2 * bn[2] * qixyz;
170 fkdx += inck.x;
171 fkdy += inck.y;
172 fkdz += inck.z;
173
174 // p-field
175 if (pscale != dscale) {
176 if CONSTEXPR (eq<ETYP, EWALD>()) {
177 bn[1] = bn1 - (1 - pscale * scale3) * rr3;
178 bn[2] = bn2 - (1 - pscale * scale5) * rr5;
179 bn[3] = bn3 - (1 - pscale * scale7) * rr7;
180 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
181 bn[1] = pscale * scale3 * rr3;
182 bn[2] = pscale * scale5 * rr5;
183 bn[3] = pscale * scale7 * rr7;
184 }
185
186 c1 = -(bn[1] * ck - bn[2] * dkr + bn[3] * qkr);
187 inci = c1 * dr - bn[1] * dkxyz + 2 * bn[2] * qkxyz;
188
189 c1 = (bn[1] * ci + bn[2] * dir + bn[3] * qir);
190 inck = c1 * dr - bn[1] * dixyz - 2 * bn[2] * qixyz;
191 }
192 fipx += inci.x;
193 fipy += inci.y;
194 fipz += inci.z;
195 fkpx += inck.x;
196 fkpy += inck.y;
197 fkpz += inck.z;
198}
199
200#pragma acc routine seq
201template <class ETYP>
203void pair_ufield(real r2, real xr, real yr, real zr, real uscale, //
204 real uindi0, real uindi1, real uindi2, real uinpi0, real uinpi1, real uinpi2,
205 real pdi,
206 real pti, //
207 real uindk0, real uindk1, real uindk2, real uinpk0, real uinpk1, real uinpk2,
208 real pdk,
209 real ptk, //
210 real aewald, real3& restrict fid, real3& restrict fip, real3& restrict fkd,
211 real3& restrict fkp)
212{
213 real r = REAL_SQRT(r2);
214 real invr1 = REAL_RECIP(r);
215 real rr2 = invr1 * invr1;
216
217 real scale3, scale5;
218 damp_thole2(r, pdi, pti, pdk, ptk, scale3, scale5);
219
220 real bn[3];
221 if CONSTEXPR (eq<ETYP, EWALD>()) damp_ewald<3>(bn, r, invr1, rr2, aewald);
222 real rr1 = invr1;
223 real rr3 = rr1 * rr2;
224 real rr5 = 3 * rr1 * rr2 * rr2;
225
226 if CONSTEXPR (eq<ETYP, EWALD>()) {
227 bn[1] -= (1 - scale3) * rr3;
228 bn[2] -= (1 - scale5) * rr5;
229 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
230 bn[1] = uscale * scale3 * rr3;
231 bn[2] = uscale * scale5 * rr5;
232 }
233
234 real coef;
235 real3 dr = make_real3(xr, yr, zr);
236 real3 uid = make_real3(uindi0, uindi1, uindi2);
237 real3 uip = make_real3(uinpi0, uinpi1, uinpi2);
238 real3 ukd = make_real3(uindk0, uindk1, uindk2);
239 real3 ukp = make_real3(uinpk0, uinpk1, uinpk2);
240
241 coef = bn[2] * dot3(dr, ukd);
242 fid += coef * dr - bn[1] * ukd;
243
244 coef = bn[2] * dot3(dr, ukp);
245 fip += coef * dr - bn[1] * ukp;
246
247 coef = bn[2] * dot3(dr, uid);
248 fkd += coef * dr - bn[1] * uid;
249
250 coef = bn[2] * dot3(dr, uip);
251 fkp += coef * dr - bn[1] * uip;
252}
253
254#pragma acc routine seq
255template <class ETYP>
257void pair_ufield_v2(real r2, real xr, real yr, real zr, real uscale,
258 real aewald, //
259 real uindi0, real uindi1, real uindi2, real uinpi0, real uinpi1, real uinpi2,
260 real pdi,
261 real pti, //
262 real uindk0, real uindk1, real uindk2, real uinpk0, real uinpk1, real uinpk2,
263 real pdk,
264 real ptk, //
265 real& restrict fidx, real& restrict fidy, real& restrict fidz,
266 real& restrict fipx, real& restrict fipy, real& restrict fipz,
267 real& restrict fkdx, real& restrict fkdy, real& restrict fkdz,
268 real& restrict fkpx, real& restrict fkpy, real& restrict fkpz)
269{
270 real r = REAL_SQRT(r2);
271 real invr1 = REAL_RECIP(r);
272 real rr2 = invr1 * invr1;
273
274 real scale3, scale5;
275 damp_thole2(r, pdi, pti, pdk, ptk, scale3, scale5);
276
277 real bn[3];
278 if CONSTEXPR (eq<ETYP, EWALD>()) damp_ewald<3>(bn, r, invr1, rr2, aewald);
279 real rr1 = invr1;
280 real rr3 = rr1 * rr2;
281 real rr5 = 3 * rr1 * rr2 * rr2;
282
283 if CONSTEXPR (eq<ETYP, EWALD>()) {
284 bn[1] -= (1 - uscale * scale3) * rr3;
285 bn[2] -= (1 - uscale * scale5) * rr5;
286 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
287 bn[1] = uscale * scale3 * rr3;
288 bn[2] = uscale * scale5 * rr5;
289 }
290
291 real coef;
292 real3 dr = make_real3(xr, yr, zr);
293 real3 uid = make_real3(uindi0, uindi1, uindi2);
294 real3 uip = make_real3(uinpi0, uinpi1, uinpi2);
295 real3 ukd = make_real3(uindk0, uindk1, uindk2);
296 real3 ukp = make_real3(uinpk0, uinpk1, uinpk2);
297
298 coef = bn[2] * dot3(dr, ukd);
299 real3 fid = coef * dr - bn[1] * ukd;
300 fidx += fid.x;
301 fidy += fid.y;
302 fidz += fid.z;
303
304 coef = bn[2] * dot3(dr, ukp);
305 real3 fip = coef * dr - bn[1] * ukp;
306 fipx += fip.x;
307 fipy += fip.y;
308 fipz += fip.z;
309
310 coef = bn[2] * dot3(dr, uid);
311 real3 fkd = coef * dr - bn[1] * uid;
312 fkdx += fkd.x;
313 fkdy += fkd.y;
314 fkdz += fkd.z;
315
316 coef = bn[2] * dot3(dr, uip);
317 real3 fkp = coef * dr - bn[1] * uip;
318 fkpx += fkp.x;
319 fkpy += fkp.y;
320 fkpz += fkp.z;
321}
322}
float x
Definition: acc/realndef.h:49
float y
Definition: acc/realndef.h:49
float z
Definition: acc/realndef.h:49
#define SEQ_CUDA
Definition: acc/seqdef.h:12
Definition: acc/realndef.h:48
#define restrict
Definition: macro.h:51
#define CONSTEXPR
Definition: macro.h:61
__device__ real dot3(real3 a, real3 b)
Definition: realn.h:84
float real
Definition: precision.h:80
Definition: testrt.h:9
__device__ void pair_ufield_v2(real r2, real xr, real yr, real zr, real uscale, real aewald, real uindi0, real uindi1, real uindi2, real uinpi0, real uinpi1, real uinpi2, real pdi, real pti, real uindk0, real uindk1, real uindk2, real uinpk0, real uinpk1, real uinpk2, real pdk, real ptk, real &__restrict__ fidx, real &__restrict__ fidy, real &__restrict__ fidz, real &__restrict__ fipx, real &__restrict__ fipy, real &__restrict__ fipz, real &__restrict__ fkdx, real &__restrict__ fkdy, real &__restrict__ fkdz, real &__restrict__ fkpx, real &__restrict__ fkpy, real &__restrict__ fkpz)
Definition: pair_field.h:257
__device__ void damp_thole3(real r, real pdi, real pti, real pdk, real ptk, real &__restrict__ scale3, real &__restrict__ scale5, real &__restrict__ scale7)
Definition: damp.h:43
__device__ void pair_dfield(real r2, real xr, real yr, real zr, real dscale, real pscale, real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real pdi, real pti, real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real pdk, real ptk, real aewald, real3 &__restrict__ fid, real3 &__restrict__ fip, real3 &__restrict__ fkd, real3 &__restrict__ fkp)
Definition: pair_field.h:11
__device__ void pair_ufield(real r2, real xr, real yr, real zr, real uscale, real uindi0, real uindi1, real uindi2, real uinpi0, real uinpi1, real uinpi2, real pdi, real pti, real uindk0, real uindk1, real uindk2, real uinpk0, real uinpk1, real uinpk2, real pdk, real ptk, real aewald, real3 &__restrict__ fid, real3 &__restrict__ fip, real3 &__restrict__ fkd, real3 &__restrict__ fkp)
Definition: pair_field.h:203
__device__ void damp_thole2(real r, real pdi, real pti, real pdk, real ptk, real &__restrict__ scale3, real &__restrict__ scale5)
Definition: damp.h:21
__device__ void pair_dfield_v2(real r2, real xr, real yr, real zr, real dscale, real pscale, real aewald, real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real pdi, real pti, real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real pdk, real ptk, real &__restrict__ fidx, real &__restrict__ fidy, real &__restrict__ fidz, real &__restrict__ fipx, real &__restrict__ fipy, real &__restrict__ fipz, real &__restrict__ fkdx, real &__restrict__ fkdy, real &__restrict__ fkdz, real &__restrict__ fkpx, real &__restrict__ fkpy, real &__restrict__ fkpz)
Definition: pair_field.h:101