Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pair_polar_chgpen.h
1#pragma once
2#include "seq/damp_hippo.h"
3#include "seq/pair_polar.h"
4
5namespace tinker {
6#pragma acc routine seq
7template <bool do_e, bool do_g, class ETYP, int CFLX>
8SEQ_CUDA
9void pair_polar_chgpen(real r2, real xr, real yr, real zr, real dscale,
10 real wscale, real ci, real dix, real diy, real diz, real corei, real vali,
11 real alphai, real qixx, real qixy, real qixz, real qiyy, real qiyz,
12 real qizz, real uix, real uiy, real uiz, real ck, real dkx, real dky,
13 real dkz, real corek, real valk, real alphak, real qkxx, real qkxy,
14 real qkxz, real qkyy, real qkyz, real qkzz, real ukx, real uky, real ukz,
15 real f, real aewald, real& restrict e, real& restrict poti,
16 real& restrict potk, PairPolarGrad& restrict pgrad)
17{
18 real dir = dix * xr + diy * yr + diz * zr;
19 real qix = qixx * xr + qixy * yr + qixz * zr;
20 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
21 real qiz = qixz * xr + qiyz * yr + qizz * zr;
22 real qir = qix * xr + qiy * yr + qiz * zr;
23 real dkr = dkx * xr + dky * yr + dkz * zr;
24 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
25 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
26 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
27 real qkr = qkx * xr + qky * yr + qkz * zr;
28 real uir = uix * xr + uiy * yr + uiz * zr;
29 real ukr = ukx * xr + uky * yr + ukz * zr;
30
31 real r = REAL_SQRT(r2);
32 real invr1 = REAL_RECIP(r);
33 real rr2 = invr1 * invr1;
34
35 real rr1 = f * invr1;
36 real rr3 = rr1 * rr2;
37 real rr5 = 3 * rr3 * rr2;
38 real rr7 = 5 * rr5 * rr2;
39 real rr9 = 7 * rr7 * rr2;
40
41 real bn[5];
42 real dmpi[5], dmpk[5];
43 real dmpik[6];
44
45 real rr3core, rr5core, rr3i, rr5i, rr7i, rr9i;
46 real rr3k, rr5k, rr7k, rr9k, rr5ik, rr7ik;
47 real dsr3i, dsr5i, dsr7i, dsr3k, dsr5k, dsr7k;
48 if CONSTEXPR (eq<ETYP, EWALD>()) {
49 if CONSTEXPR (do_g) {
50 damp_ewald<5>(bn, r, invr1, rr2, aewald);
51 damp_pole_v2<11>(dmpik, dmpi, dmpk, r, alphai, alphak);
52 } else {
53 damp_ewald<4>(bn, r, invr1, rr2, aewald);
54 damp_pole_v2<9>(dmpik, dmpi, dmpk, r, alphai, alphak);
55 }
56
57 bn[1] *= f;
58 bn[2] *= f;
59 bn[3] *= f;
60 bn[4] *= f;
61 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
62 if CONSTEXPR (do_g)
63 damp_pole_v2<11>(dmpik, dmpi, dmpk, r, alphai, alphak);
64 else {
65 damp_pole_v2<9>(dmpik, dmpi, dmpk, r, alphai, alphak);
66 }
67
68 bn[1] = rr3;
69 bn[2] = rr5;
70 bn[3] = rr7;
71 bn[4] = rr9;
72 }
73
74 rr3core = bn[1] - (1 - dscale) * rr3;
75 rr5core = bn[2] - (1 - dscale) * rr5;
76 rr3i = bn[1] - (1 - dscale * dmpi[1]) * rr3;
77 rr5i = bn[2] - (1 - dscale * dmpi[2]) * rr5;
78 rr7i = bn[3] - (1 - dscale * dmpi[3]) * rr7;
79 rr9i = bn[4] - (1 - dscale * dmpi[4]) * rr9;
80 rr3k = bn[1] - (1 - dscale * dmpk[1]) * rr3;
81 rr5k = bn[2] - (1 - dscale * dmpk[2]) * rr5;
82 rr7k = bn[3] - (1 - dscale * dmpk[3]) * rr7;
83 rr9k = bn[4] - (1 - dscale * dmpk[4]) * rr9;
84 rr5ik = bn[2] - (1 - wscale * dmpik[2]) * rr5;
85 rr7ik = bn[3] - (1 - wscale * dmpik[3]) * rr7;
86
87 dsr3i = 2 * rr3i;
88 dsr3k = 2 * rr3k;
89 dsr5i = 2 * rr5i;
90 dsr5k = 2 * rr5k;
91 dsr7i = 2 * rr7i;
92 dsr7k = 2 * rr7k;
93
94 if CONSTEXPR (do_e) {
95 real diu = dix * ukx + diy * uky + diz * ukz;
96 real qiu = qix * ukx + qiy * uky + qiz * ukz;
97 real dku = dkx * uix + dky * uiy + dkz * uiz;
98 real qku = qkx * uix + qky * uiy + qkz * uiz;
99 e = uir * (corek * rr3core + valk * rr3k)
100 - ukr * (corei * rr3core + vali * rr3i) + diu * rr3i + dku * rr3k
101 + 2 * (qiu * rr5i - qku * rr5k) - dkr * uir * rr5k - dir * ukr * rr5i
102 + qkr * uir * rr7k - qir * ukr * rr7i;
103 }
104
105 if CONSTEXPR (do_g) {
106 // get the induced dipole field used for dipole torques
107 real tix3 = dsr3i * ukx;
108 real tiy3 = dsr3i * uky;
109 real tiz3 = dsr3i * ukz;
110 real tkx3 = dsr3k * uix;
111 real tky3 = dsr3k * uiy;
112 real tkz3 = dsr3k * uiz;
113 real tuir = -dsr5i * ukr;
114 real tukr = -dsr5k * uir;
115
116 if CONSTEXPR (CFLX) {
117 poti = -ukr * dsr3i;
118 potk = uir * dsr3k;
119 }
120
121 pgrad.ufldi[0] = (tix3 + xr * tuir);
122 pgrad.ufldi[1] = (tiy3 + yr * tuir);
123 pgrad.ufldi[2] = (tiz3 + zr * tuir);
124 pgrad.ufldk[0] = (tkx3 + xr * tukr);
125 pgrad.ufldk[1] = (tky3 + yr * tukr);
126 pgrad.ufldk[2] = (tkz3 + zr * tukr);
127
128 // get induced dipole field gradient used for quadrupole torques
129
130 real tix5 = 2 * (dsr5i * ukx);
131 real tiy5 = 2 * (dsr5i * uky);
132 real tiz5 = 2 * (dsr5i * ukz);
133 real tkx5 = 2 * (dsr5k * uix);
134 real tky5 = 2 * (dsr5k * uiy);
135 real tkz5 = 2 * (dsr5k * uiz);
136 tuir = -dsr7i * ukr;
137 tukr = -dsr7k * uir;
138
139 pgrad.dufldi[0] = (xr * tix5 + xr * xr * tuir);
140 pgrad.dufldi[1] = (xr * tiy5 + yr * tix5 + 2 * xr * yr * tuir);
141 pgrad.dufldi[2] = (yr * tiy5 + yr * yr * tuir);
142 pgrad.dufldi[3] = (xr * tiz5 + zr * tix5 + 2 * xr * zr * tuir);
143 pgrad.dufldi[4] = (yr * tiz5 + zr * tiy5 + 2 * yr * zr * tuir);
144 pgrad.dufldi[5] = (zr * tiz5 + zr * zr * tuir);
145 pgrad.dufldk[0] = (-xr * tkx5 - xr * xr * tukr);
146 pgrad.dufldk[1] = (-xr * tky5 - yr * tkx5 - 2 * xr * yr * tukr);
147 pgrad.dufldk[2] = (-yr * tky5 - yr * yr * tukr);
148 pgrad.dufldk[3] = (-xr * tkz5 - zr * tkx5 - 2 * xr * zr * tukr);
149 pgrad.dufldk[4] = (-yr * tkz5 - zr * tky5 - 2 * yr * zr * tukr);
150 pgrad.dufldk[5] = (-zr * tkz5 - zr * zr * tukr);
151
152 // get the field gradient for direct polarization force
153
154 real term1i, term2i, term3i, term4i, term5i, term6i, term1core;
155 real term7i, term8i;
156 real term1k, term2k, term3k, term4k, term5k, term6k;
157 real term7k, term8k;
158
159 term1i = rr3i - rr5i * xr * xr;
160 term1core = rr3core - rr5core * xr * xr;
161 term2i = 2 * rr5i * xr;
162 term3i = rr7i * xr * xr - rr5i;
163 term4i = 2 * rr5i;
164 term5i = 5 * rr7i * xr;
165 term6i = rr9i * xr * xr;
166 term1k = rr3k - rr5k * xr * xr;
167 term2k = 2 * rr5k * xr;
168 term3k = rr7k * xr * xr - rr5k;
169 term4k = 2 * rr5k;
170 term5k = 5 * rr7k * xr;
171 term6k = rr9k * xr * xr;
172 real tixx = vali * term1i + corei * term1core + dix * term2i
173 - dir * term3i - qixx * term4i + qix * term5i - qir * term6i
174 + (qiy * yr + qiz * zr) * rr7i;
175 real tkxx = valk * term1k + corek * term1core - dkx * term2k
176 + dkr * term3k - qkxx * term4k + qkx * term5k - qkr * term6k
177 + (qky * yr + qkz * zr) * rr7k;
178 term1i = rr3i - rr5i * yr * yr;
179 term1core = rr3core - rr5core * yr * yr;
180 term2i = 2 * rr5i * yr;
181 term3i = rr7i * yr * yr - rr5i;
182 term4i = 2 * rr5i;
183 term5i = 5 * rr7i * yr;
184 term6i = rr9i * yr * yr;
185 term1k = rr3k - rr5k * yr * yr;
186 term2k = 2 * rr5k * yr;
187 term3k = rr7k * yr * yr - rr5k;
188 term4k = 2 * rr5k;
189 term5k = 5 * rr7k * yr;
190 term6k = rr9k * yr * yr;
191 real tiyy = vali * term1i + corei * term1core + diy * term2i
192 - dir * term3i - qiyy * term4i + qiy * term5i - qir * term6i
193 + (qix * xr + qiz * zr) * rr7i;
194 real tkyy = valk * term1k + corek * term1core - dky * term2k
195 + dkr * term3k - qkyy * term4k + qky * term5k - qkr * term6k
196 + (qkx * xr + qkz * zr) * rr7k;
197 term1i = rr3i - rr5i * zr * zr;
198 term1core = rr3core - rr5core * zr * zr;
199 term2i = 2 * rr5i * zr;
200 term3i = rr7i * zr * zr - rr5i;
201 term4i = 2 * rr5i;
202 term5i = 5 * rr7i * zr;
203 term6i = rr9i * zr * zr;
204 term1k = rr3k - rr5k * zr * zr;
205 term2k = 2 * rr5k * zr;
206 term3k = rr7k * zr * zr - rr5k;
207 term4k = 2 * rr5k;
208 term5k = 5 * rr7k * zr;
209 term6k = rr9k * zr * zr;
210 real tizz = vali * term1i + corei * term1core + diz * term2i
211 - dir * term3i - qizz * term4i + qiz * term5i - qir * term6i
212 + (qix * xr + qiy * yr) * rr7i;
213 real tkzz = valk * term1k + corek * term1core - dkz * term2k
214 + dkr * term3k - qkzz * term4k + qkz * term5k - qkr * term6k
215 + (qkx * xr + qky * yr) * rr7k;
216 term2i = rr5i * xr;
217 term1i = yr * term2i;
218 term1core = rr5core * xr * yr;
219 term3i = rr5i * yr;
220 term4i = yr * (rr7i * xr);
221 term5i = 2 * rr5i;
222 term6i = 2 * rr7i * xr;
223 term7i = 2 * rr7i * yr;
224 term8i = yr * rr9i * xr;
225 term2k = rr5k * xr;
226 term1k = yr * term2k;
227 term3k = rr5k * yr;
228 term4k = yr * (rr7k * xr);
229 term5k = 2 * rr5k;
230 term6k = 2 * rr7k * xr;
231 term7k = 2 * rr7k * yr;
232 term8k = yr * rr9k * xr;
233 real tixy = -vali * term1i - corei * term1core + diy * term2i
234 + dix * term3i - dir * term4i - qixy * term5i + qiy * term6i
235 + qix * term7i - qir * term8i;
236 real tkxy = -valk * term1k - corek * term1core - dky * term2k
237 - dkx * term3k + dkr * term4k - qkxy * term5k + qky * term6k
238 + qkx * term7k - qkr * term8k;
239
240 term2i = rr5i * xr;
241 term1i = zr * term2i;
242 term1core = rr5core * xr * zr;
243 term3i = rr5i * zr;
244 term4i = zr * (rr7i * xr);
245 term5i = 2 * rr5i;
246 term6i = 2 * rr7i * xr;
247 term7i = 2 * rr7i * zr;
248 term8i = zr * rr9i * xr;
249 term2k = rr5k * xr;
250 term1k = zr * term2k;
251 term3k = rr5k * zr;
252 term4k = zr * (rr7k * xr);
253 term5k = 2 * rr5k;
254 term6k = 2 * rr7k * xr;
255 term7k = 2 * rr7k * zr;
256 term8k = zr * rr9k * xr;
257 real tixz = -vali * term1i - corei * term1core + diz * term2i
258 + dix * term3i - dir * term4i - qixz * term5i + qiz * term6i
259 + qix * term7i - qir * term8i;
260 real tkxz = -valk * term1k - corek * term1core - dkz * term2k
261 - dkx * term3k + dkr * term4k - qkxz * term5k + qkz * term6k
262 + qkx * term7k - qkr * term8k;
263
264 term2i = rr5i * yr;
265 term1i = zr * term2i;
266 term1core = rr5core * yr * zr;
267 term3i = rr5i * zr;
268 term4i = zr * (rr7i * yr);
269 term5i = 2 * rr5i;
270 term6i = 2 * rr7i * yr;
271 term7i = 2 * rr7i * zr;
272 term8i = zr * rr9i * yr;
273 term2k = rr5k * yr;
274 term1k = zr * term2k;
275 term3k = rr5k * zr;
276 term4k = zr * (rr7k * yr);
277 term5k = 2 * rr5k;
278 term6k = 2 * rr7k * yr;
279 term7k = 2 * rr7k * zr;
280 term8k = zr * rr9k * yr;
281 real tiyz = -vali * term1i - corei * term1core + diz * term2i
282 + diy * term3i - dir * term4i - qiyz * term5i + qiz * term6i
283 + qiy * term7i - qir * term8i;
284 real tkyz = -valk * term1k - corek * term1core - dkz * term2k
285 - dky * term3k + dkr * term4k - qkyz * term5k + qkz * term6k
286 + qky * term7k - qkr * term8k;
287
288 real depx, depy, depz;
289
290 depx = tixx * ukx + tixy * uky + tixz * ukz - tkxx * uix - tkxy * uiy
291 - tkxz * uiz;
292 depy = tixy * ukx + tiyy * uky + tiyz * ukz - tkxy * uix - tkyy * uiy
293 - tkyz * uiz;
294 depz = tixz * ukx + tiyz * uky + tizz * ukz - tkxz * uix - tkyz * uiy
295 - tkzz * uiz;
296
297 pgrad.frcx = 2 * depx;
298 pgrad.frcy = 2 * depy;
299 pgrad.frcz = 2 * depz;
300
301 // get the dtau/dr terms used for mutual polarization force
302
303 real term1 = 2 * rr5ik;
304 real term2 = term1 * xr;
305 real term3 = rr5ik - rr7ik * xr * xr;
306 tixx = uix * term2 + uir * term3;
307 tkxx = ukx * term2 + ukr * term3;
308 term2 = term1 * yr;
309 term3 = rr5ik - rr7ik * yr * yr;
310 tiyy = uiy * term2 + uir * term3;
311 tkyy = uky * term2 + ukr * term3;
312 term2 = term1 * zr;
313 term3 = rr5ik - rr7ik * zr * zr;
314 tizz = uiz * term2 + uir * term3;
315 tkzz = ukz * term2 + ukr * term3;
316 term1 = rr5ik * yr;
317 term2 = rr5ik * xr;
318 term3 = yr * (rr7ik * xr);
319 tixy = uix * term1 + uiy * term2 - uir * term3;
320 tkxy = ukx * term1 + uky * term2 - ukr * term3;
321 term1 = rr5ik * zr;
322 term3 = zr * (rr7ik * xr);
323 tixz = uix * term1 + uiz * term2 - uir * term3;
324 tkxz = ukx * term1 + ukz * term2 - ukr * term3;
325 term2 = rr5ik * yr;
326 term3 = zr * (rr7ik * yr);
327 tiyz = uiy * term1 + uiz * term2 - uir * term3;
328 tkyz = uky * term1 + ukz * term2 - ukr * term3;
329
330 depx = tixx * ukx + tixy * uky + tixz * ukz + tkxx * uix + tkxy * uiy
331 + tkxz * uiz;
332 depy = tixy * ukx + tiyy * uky + tiyz * ukz + tkxy * uix + tkyy * uiy
333 + tkyz * uiz;
334 depz = tixz * ukx + tiyz * uky + tizz * ukz + tkxz * uix + tkyz * uiy
335 + tkzz * uiz;
336
337 pgrad.frcx += depx;
338 pgrad.frcy += depy;
339 pgrad.frcz += depz;
340 }
341}
342}
SEQ_CUDA
#define SEQ_CUDA
Definition: acc/seqdef.h:12
restrict
#define restrict
Definition: macro.h:51
CONSTEXPR
#define CONSTEXPR
Definition: macro.h:61
tinker::real
float real
Definition: precision.h:80
tinker
Definition: testrt.h:9
tinker::depy
grad_prec * depy
tinker::pair_polar_chgpen
__device__ void pair_polar_chgpen(real r2, real xr, real yr, real zr, real dscale, real wscale, real ci, real dix, real diy, real diz, real corei, real vali, real alphai, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real uix, real uiy, real uiz, real ck, real dkx, real dky, real dkz, real corek, real valk, real alphak, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real ukx, real uky, real ukz, real f, real aewald, real &__restrict__ e, real &__restrict__ poti, real &__restrict__ potk, PairPolarGrad &__restrict__ pgrad)
Definition: pair_polar_chgpen.h:9
tinker::depx
grad_prec * depx
tinker::depz
grad_prec * depz
tinker::PairPolarGrad
Definition: pair_polar.h:8
Generated by doxygen 1.9.5