Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pairpolaraplus.h
1#pragma once
2#include "seq/damp.h"
3#include "seq/dampaplus.h"
4#include "seq/pair_polar.h"
5
6namespace tinker {
7#pragma acc routine seq
8template <bool do_e, bool do_g, class ETYP, int CFLX>
9SEQ_CUDA
10void pair_polar_aplus(real r2, real xr, real yr, real zr, real dscale,
11 real uscale, real ci, real dix, real diy, real diz, real pdi, real pti,
12 real ddi, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz,
13 real uix, real uiy, real uiz, real ck, real dkx, real dky, real dkz,
14 real pdk, real ptk, real ddk, real qkxx, real qkxy, real qkxz, real qkyy,
15 real qkyz, real qkzz, real ukx, real uky, real ukz, real f, real aewald,
16 real& restrict e, real& restrict poti, real& restrict potk,
17 PairPolarGrad& restrict pgrad)
18{
19 real dir = dix * xr + diy * yr + diz * zr;
20 real qix = qixx * xr + qixy * yr + qixz * zr;
21 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
22 real qiz = qixz * xr + qiyz * yr + qizz * zr;
23 real qir = qix * xr + qiy * yr + qiz * zr;
24 real dkr = dkx * xr + dky * yr + dkz * zr;
25 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
26 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
27 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
28 real qkr = qkx * xr + qky * yr + qkz * zr;
29 real uir = uix * xr + uiy * yr + uiz * zr;
30 real ukr = ukx * xr + uky * yr + ukz * zr;
31
32 real r = REAL_SQRT(r2);
33 real invr1 = REAL_RECIP(r);
34 real rr2 = invr1 * invr1;
35
36 real rr1 = invr1;
37 real rr3 = rr1 * rr2;
38 real rr5 = 3 * rr3 * rr2;
39 real rr7 = 5 * rr5 * rr2;
40 MAYBE_UNUSED real rr9;
41 if CONSTEXPR (do_g) rr9 = 7 * rr7 * rr2;
42 real bn[5];
43 if CONSTEXPR (eq<ETYP, EWALD>()) {
44 if CONSTEXPR (!do_g)
45 damp_ewald<4>(bn, r, invr1, rr2, aewald);
46 else
47 damp_ewald<5>(bn, r, invr1, rr2, aewald);
48 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
49 bn[1] = rr3;
50 bn[2] = rr5;
51 bn[3] = rr7;
52 if CONSTEXPR (do_g) bn[4] = rr9;
53 }
54
55 // if use_dirdamp
56 real ex3, ex5, ex7;
57 MAYBE_UNUSED real rc31, rc32, rc33, rc51, rc52, rc53, rc71, rc72, rc73;
58 if CONSTEXPR (!do_g) {
59 damp_aplus3(r, pdi, ddi, pdk, ddk, //
60 ex3, ex5, ex7);
61 ex3 = 1 - ex3;
62 ex5 = 1 - ex5;
63 ex7 = 1 - ex7;
64 } else {
65 damp_aplus3g( //
66 r, rr2, xr, yr, zr, //
67 pdi, ddi, pdk, ddk, //
68 ex3, ex5, ex7, //
69 rc31, rc32, rc33, //
70 rc51, rc52, rc53, //
71 rc71, rc72, rc73);
72 rc31 *= rr3;
73 rc32 *= rr3;
74 rc33 *= rr3;
75 rc51 *= rr5;
76 rc52 *= rr5;
77 rc53 *= rr5;
78 rc71 *= rr7;
79 rc72 *= rr7;
80 rc73 *= rr7;
81 }
82 // end if use_dirdamp
83
84 real sr3 = bn[1] - ex3 * rr3;
85 real sr5 = bn[2] - ex5 * rr5;
86 real sr7 = bn[3] - ex7 * rr7;
87
88 if CONSTEXPR (do_e) {
89 real diu = dix * ukx + diy * uky + diz * ukz;
90 real qiu = qix * ukx + qiy * uky + qiz * ukz;
91 real dku = dkx * uix + dky * uiy + dkz * uiz;
92 real qku = qkx * uix + qky * uiy + qkz * uiz;
93 real term1 = ck * uir - ci * ukr + diu + dku;
94 real term2 = 2 * (qiu - qku) - uir * dkr - dir * ukr;
95 real term3 = uir * qkr - ukr * qir;
96 e = dscale * f * (term1 * sr3 + term2 * sr5 + term3 * sr7);
97 }
98
99 if CONSTEXPR (do_g) {
100 real uir = uix * xr + uiy * yr + uiz * zr;
101 real ukr = ukx * xr + uky * yr + ukz * zr;
102 if CONSTEXPR (CFLX) {
103 poti = -2.0 * ukr * f * dscale * sr3;
104 potk = 2.0 * uir * f * dscale * sr3;
105 }
106
107 // get the induced dipole field used for dipole torques
108
109 real tuir, tukr;
110
111 real tix3 = 2.0 * dscale * sr3 * ukx;
112 real tiy3 = 2.0 * dscale * sr3 * uky;
113 real tiz3 = 2.0 * dscale * sr3 * ukz;
114 real tkx3 = 2.0 * dscale * sr3 * uix;
115 real tky3 = 2.0 * dscale * sr3 * uiy;
116 real tkz3 = 2.0 * dscale * sr3 * uiz;
117 tuir = -2.0 * dscale * sr5 * ukr;
118 tukr = -2.0 * dscale * sr5 * uir;
119
120 pgrad.ufldi[0] = f * (tix3 + xr * tuir);
121 pgrad.ufldi[1] = f * (tiy3 + yr * tuir);
122 pgrad.ufldi[2] = f * (tiz3 + zr * tuir);
123 pgrad.ufldk[0] = f * (tkx3 + xr * tukr);
124 pgrad.ufldk[1] = f * (tky3 + yr * tukr);
125 pgrad.ufldk[2] = f * (tkz3 + zr * tukr);
126
127 // get induced dipole field gradient used for quadrupole torques
128
129 real tix5 = 4 * (dscale * sr5 * ukx);
130 real tiy5 = 4 * (dscale * sr5 * uky);
131 real tiz5 = 4 * (dscale * sr5 * ukz);
132 real tkx5 = 4 * (dscale * sr5 * uix);
133 real tky5 = 4 * (dscale * sr5 * uiy);
134 real tkz5 = 4 * (dscale * sr5 * uiz);
135 tuir = -2.0 * dscale * sr7 * ukr;
136 tukr = -2.0 * dscale * sr7 * uir;
137
138 pgrad.dufldi[0] = f * (xr * tix5 + xr * xr * tuir);
139 pgrad.dufldi[1] = f * (xr * tiy5 + yr * tix5 + 2 * xr * yr * tuir);
140 pgrad.dufldi[2] = f * (yr * tiy5 + yr * yr * tuir);
141 pgrad.dufldi[3] = f * (xr * tiz5 + zr * tix5 + 2 * xr * zr * tuir);
142 pgrad.dufldi[4] = f * (yr * tiz5 + zr * tiy5 + 2 * yr * zr * tuir);
143 pgrad.dufldi[5] = f * (zr * tiz5 + zr * zr * tuir);
144 pgrad.dufldk[0] = f * (-xr * tkx5 - xr * xr * tukr);
145 pgrad.dufldk[1] = f * (-xr * tky5 - yr * tkx5 - 2 * xr * yr * tukr);
146 pgrad.dufldk[2] = f * (-yr * tky5 - yr * yr * tukr);
147 pgrad.dufldk[3] = f * (-xr * tkz5 - zr * tkx5 - 2 * xr * zr * tukr);
148 pgrad.dufldk[4] = f * (-yr * tkz5 - zr * tky5 - 2 * yr * zr * tukr);
149 pgrad.dufldk[5] = f * (-zr * tkz5 - zr * zr * tukr);
150
151 // get the field gradient for direct polarization force
152
153 real term1, term2, term3, term4, term5, term6, term7;
154
155 term1 = bn[2] - ex3 * rr5;
156 term2 = bn[3] - ex5 * rr7;
157 term3 = -sr3 + term1 * xr * xr - xr * rc31;
158 term4 = rc31 - term1 * xr - sr5 * xr;
159 term5 = term2 * xr * xr - sr5 - xr * rc51;
160 term6 = (bn[4] - ex7 * rr9) * xr * xr - bn[3] - xr * rc71;
161 term7 = rc51 - 2 * bn[3] * xr + (ex5 + 1.5f * ex7) * rr7 * xr;
162 real tixx = ci * term3 + dix * term4 + dir * term5 + 2 * sr5 * qixx
163 + (qiy * yr + qiz * zr) * ex7 * rr7 + 2 * qix * term7 + qir * term6;
164 real tkxx = ck * term3 - dkx * term4 - dkr * term5 + 2 * sr5 * qkxx
165 + (qky * yr + qkz * zr) * ex7 * rr7 + 2 * qkx * term7 + qkr * term6;
166
167 term3 = -sr3 + term1 * yr * yr - yr * rc32;
168 term4 = rc32 - term1 * yr - sr5 * yr;
169 term5 = term2 * yr * yr - sr5 - yr * rc52;
170 term6 = (bn[4] - ex7 * rr9) * yr * yr - bn[3] - yr * rc72;
171 term7 = rc52 - 2 * bn[3] * yr + (ex5 + 1.5f * ex7) * rr7 * yr;
172 real tiyy = ci * term3 + diy * term4 + dir * term5 + 2 * sr5 * qiyy
173 + (qix * xr + qiz * zr) * ex7 * rr7 + 2 * qiy * term7 + qir * term6;
174 real tkyy = ck * term3 - dky * term4 - dkr * term5 + 2 * sr5 * qkyy
175 + (qkx * xr + qkz * zr) * ex7 * rr7 + 2 * qky * term7 + qkr * term6;
176
177 term3 = -sr3 + term1 * zr * zr - zr * rc33;
178 term4 = rc33 - term1 * zr - sr5 * zr;
179 term5 = term2 * zr * zr - sr5 - zr * rc53;
180 term6 = (bn[4] - ex7 * rr9) * zr * zr - bn[3] - zr * rc73;
181 term7 = rc53 - 2 * bn[3] * zr + (ex5 + 1.5f * ex7) * rr7 * zr;
182 real tizz = ci * term3 + diz * term4 + dir * term5 + 2 * sr5 * qizz
183 + (qix * xr + qiy * yr) * ex7 * rr7 + 2 * qiz * term7 + qir * term6;
184 real tkzz = ck * term3 - dkz * term4 - dkr * term5 + 2 * sr5 * qkzz
185 + (qkx * xr + qky * yr) * ex7 * rr7 + 2 * qkz * term7 + qkr * term6;
186
187 term3 = term1 * xr * yr - yr * rc31;
188 term4 = rc31 - term1 * xr;
189 term5 = term2 * xr * yr - yr * rc51;
190 term6 = (bn[4] - ex7 * rr9) * xr * yr - yr * rc71;
191 term7 = rc51 - term2 * xr;
192 real tixy = ci * term3 - sr5 * dix * yr + diy * term4 + dir * term5
193 + 2 * sr5 * qixy - 2 * sr7 * yr * qix + 2 * qiy * term7 + qir * term6;
194 real tkxy = ck * term3 + sr5 * dkx * yr - dky * term4 - dkr * term5
195 + 2 * sr5 * qkxy - 2 * sr7 * yr * qkx + 2 * qky * term7 + qkr * term6;
196
197 term3 = term1 * xr * zr - zr * rc31;
198 term5 = term2 * xr * zr - zr * rc51;
199 term6 = (bn[4] - ex7 * rr9) * xr * zr - zr * rc71;
200 real tixz = ci * term3 - sr5 * dix * zr + diz * term4 + dir * term5
201 + 2 * sr5 * qixz - 2 * sr7 * zr * qix + 2 * qiz * term7 + qir * term6;
202 real tkxz = ck * term3 + sr5 * dkx * zr - dkz * term4 - dkr * term5
203 + 2 * sr5 * qkxz - 2 * sr7 * zr * qkx + 2 * qkz * term7 + qkr * term6;
204
205 term3 = term1 * yr * zr - zr * rc32;
206 term4 = rc32 - term1 * yr;
207 term5 = term2 * yr * zr - zr * rc52;
208 term6 = (bn[4] - ex7 * rr9) * yr * zr - zr * rc72;
209 term7 = rc52 - term2 * yr;
210 real tiyz = ci * term3 - sr5 * diy * zr + diz * term4 + dir * term5
211 + 2 * sr5 * qiyz - 2 * sr7 * zr * qiy + 2 * qiz * term7 + qir * term6;
212 real tkyz = ck * term3 + sr5 * dky * zr - dkz * term4 - dkr * term5
213 + 2 * sr5 * qkyz - 2 * sr7 * zr * qky + 2 * qkz * term7 + qkr * term6;
214
215 // get the dEd/dR terms for Thole direct polarization force
216
217 real depx, depy, depz;
218
219 depx = tixx * ukx + tixy * uky + tixz * ukz - tkxx * uix - tkxy * uiy
220 - tkxz * uiz;
221 depy = tixy * ukx + tiyy * uky + tiyz * ukz - tkxy * uix - tkyy * uiy
222 - tkyz * uiz;
223 depz = tixz * ukx + tiyz * uky + tizz * ukz - tkxz * uix - tkyz * uiy
224 - tkzz * uiz;
225 if CONSTEXPR (eq<ETYP, EWALD>()) {
226 pgrad.frcx = -2 * depx;
227 pgrad.frcy = -2 * depy;
228 pgrad.frcz = -2 * depz;
229 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
230 pgrad.frcx = -2 * depx * dscale;
231 pgrad.frcy = -2 * depy * dscale;
232 pgrad.frcz = -2 * depz * dscale;
233 }
234
235 // get the dtau/dr terms used for mutual polarization force
236 damp_thole3g( //
237 r, rr2, xr, yr, zr, //
238 pdi, pti, pdk, ptk, //
239 ex3, ex5, ex7, //
240 rc31, rc32, rc33, //
241 rc51, rc52, rc53, //
242 rc71, rc72, rc73);
243 rc31 *= rr3;
244 rc32 *= rr3;
245 rc33 *= rr3;
246 rc51 *= rr5;
247 rc52 *= rr5;
248 rc53 *= rr5;
249 sr3 = bn[1] - ex3 * rr3;
250 sr5 = bn[2] - ex5 * rr5;
251
252 term1 = bn[2] - ex3 * rr5;
253 term2 = bn[3] - ex5 * rr7;
254 term3 = sr5 + term1;
255
256 term5 = -xr * term3 + rc31;
257 term6 = -sr5 + xr * xr * term2 - xr * rc51;
258 tixx = uix * term5 + uir * term6;
259 tkxx = ukx * term5 + ukr * term6;
260
261 term5 = -yr * term3 + rc32;
262 term6 = -sr5 + yr * yr * term2 - yr * rc52;
263 tiyy = uiy * term5 + uir * term6;
264 tkyy = uky * term5 + ukr * term6;
265
266 term5 = -zr * term3 + rc33;
267 term6 = -sr5 + zr * zr * term2 - zr * rc53;
268 tizz = uiz * term5 + uir * term6;
269 tkzz = ukz * term5 + ukr * term6;
270
271 term4 = -sr5 * yr;
272 term5 = -xr * term1 + rc31;
273 term6 = xr * yr * term2 - yr * rc51;
274 tixy = uix * term4 + uiy * term5 + uir * term6;
275 tkxy = ukx * term4 + uky * term5 + ukr * term6;
276
277 term4 = -sr5 * zr;
278 term6 = xr * zr * term2 - zr * rc51;
279 tixz = uix * term4 + uiz * term5 + uir * term6;
280 tkxz = ukx * term4 + ukz * term5 + ukr * term6;
281
282 term5 = -yr * term1 + rc32;
283 term6 = yr * zr * term2 - zr * rc52;
284 tiyz = uiy * term4 + uiz * term5 + uir * term6;
285 tkyz = uky * term4 + ukz * term5 + ukr * term6;
286
287 depx = tixx * ukx + tixy * uky + tixz * ukz + tkxx * uix + tkxy * uiy
288 + tkxz * uiz;
289 depy = tixy * ukx + tiyy * uky + tiyz * ukz + tkxy * uix + tkyy * uiy
290 + tkyz * uiz;
291 depz = tixz * ukx + tiyz * uky + tizz * ukz + tkxz * uix + tkyz * uiy
292 + tkzz * uiz;
293 if CONSTEXPR (eq<ETYP, EWALD>()) {
294 pgrad.frcx -= depx;
295 pgrad.frcy -= depy;
296 pgrad.frcz -= depz;
297 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
298 pgrad.frcx -= uscale * depx;
299 pgrad.frcy -= uscale * depy;
300 pgrad.frcz -= uscale * depz;
301 }
302
303 pgrad.frcx *= f;
304 pgrad.frcy *= f;
305 pgrad.frcz *= f;
306 }
307}
308
309#pragma acc routine seq
310template <class Ver, class ETYP, int CFLX>
311SEQ_CUDA
312void pair_polar_aplus_v2(real r2, real xr, real yr, real zr, real dscale,
313 real uscale, real ci, real dix, real diy, real diz, real qixx, real qixy,
314 real qixz, real qiyy, real qiyz, real qizz, real uix, real uiy, real uiz,
315 real pdi, real pti, real ddi, real ck, real dkx, real dky, real dkz,
316 real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real ukx,
317 real uky, real ukz, real pdk, real ptk, real ddk, real f, real aewald,
318 real& restrict frcxi, real& restrict frcyi, real& restrict frczi,
319 real& restrict frcxk, real& restrict frcyk,
320 real& restrict frczk, //
321 real& restrict uf0i, real& restrict uf1i, real& restrict uf2i,
322 real& restrict uf0k, real& restrict uf1k, real& restrict uf2k,
323 real& restrict duf0i, real& restrict duf1i, real& restrict duf2i,
324 real& restrict duf3i, real& restrict duf4i, real& restrict duf5i,
325 real& restrict duf0k, real& restrict duf1k, real& restrict duf2k,
326 real& restrict duf3k, real& restrict duf4k, real& restrict duf5k,
327 real& restrict e, real& restrict vxx, real& restrict vxy, real& restrict vxz,
328 real& restrict vyy, real& restrict vyz, real& restrict vzz,
329 real& restrict poti, real& restrict potk)
330{
331 constexpr bool do_e = Ver::e;
332 constexpr bool do_g = Ver::g;
333 constexpr bool do_v = Ver::v;
334
335 if CONSTEXPR (eq<ETYP, EWALD>()) {
336 dscale = 1;
337 // pscale = 1;
338 uscale = 1;
339 }
340
341 real dir = dix * xr + diy * yr + diz * zr;
342 real qix = qixx * xr + qixy * yr + qixz * zr;
343 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
344 real qiz = qixz * xr + qiyz * yr + qizz * zr;
345 real qir = qix * xr + qiy * yr + qiz * zr;
346 real dkr = dkx * xr + dky * yr + dkz * zr;
347 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
348 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
349 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
350 real qkr = qkx * xr + qky * yr + qkz * zr;
351 real uir = uix * xr + uiy * yr + uiz * zr;
352 real ukr = ukx * xr + uky * yr + ukz * zr;
353
354 real r = REAL_SQRT(r2);
355 real invr1 = REAL_RECIP(r);
356 real rr2 = invr1 * invr1;
357
358 real rr1 = invr1;
359 real rr3 = rr1 * rr2;
360 real rr5 = 3 * rr3 * rr2;
361 real rr7 = 5 * rr5 * rr2;
362 MAYBE_UNUSED real rr9;
363 if CONSTEXPR (do_g) rr9 = 7 * rr7 * rr2;
364 real bn[5];
365 if CONSTEXPR (eq<ETYP, EWALD>()) {
366 if CONSTEXPR (!do_g)
367 damp_ewald<4>(bn, r, invr1, rr2, aewald);
368 else
369 damp_ewald<5>(bn, r, invr1, rr2, aewald);
370 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
371 bn[1] = rr3;
372 bn[2] = rr5;
373 bn[3] = rr7;
374 if CONSTEXPR (do_g) bn[4] = rr9;
375 }
376
377 // if use_dirdamp
378 real ex3, ex5, ex7;
379 MAYBE_UNUSED real rc31, rc32, rc33, rc51, rc52, rc53, rc71, rc72, rc73;
380 if CONSTEXPR (!do_g) {
381 damp_aplus3(r, pdi, ddi, pdk, ddk, //
382 ex3, ex5, ex7);
383 ex3 = 1 - ex3;
384 ex5 = 1 - ex5;
385 ex7 = 1 - ex7;
386 } else {
387 damp_aplus3g( //
388 r, rr2, xr, yr, zr, //
389 pdi, ddi, pdk, ddk, //
390 ex3, ex5, ex7, //
391 rc31, rc32, rc33, //
392 rc51, rc52, rc53, //
393 rc71, rc72, rc73);
394 rc31 *= rr3;
395 rc32 *= rr3;
396 rc33 *= rr3;
397 rc51 *= rr5;
398 rc52 *= rr5;
399 rc53 *= rr5;
400 rc71 *= rr7;
401 rc72 *= rr7;
402 rc73 *= rr7;
403 }
404 // end if use_dirdamp
405
406 real sr3 = bn[1] - ex3 * rr3;
407 real sr5 = bn[2] - ex5 * rr5;
408 real sr7 = bn[3] - ex7 * rr7;
409
410 if CONSTEXPR (do_e) {
411 real diu = dix * ukx + diy * uky + diz * ukz;
412 real qiu = qix * ukx + qiy * uky + qiz * ukz;
413 real dku = dkx * uix + dky * uiy + dkz * uiz;
414 real qku = qkx * uix + qky * uiy + qkz * uiz;
415 real term1 = ck * uir - ci * ukr + diu + dku;
416 real term2 = 2 * (qiu - qku) - uir * dkr - dir * ukr;
417 real term3 = uir * qkr - ukr * qir;
418 e = dscale * f * (term1 * sr3 + term2 * sr5 + term3 * sr7);
419 }
420
421 if CONSTEXPR (CFLX) {
422 poti = -2.0 * ukr * f * dscale * sr3;
423 potk = 2.0 * uir * f * dscale * sr3;
424 }
425
426 if CONSTEXPR (do_g) {
427 real uir = uix * xr + uiy * yr + uiz * zr;
428 real ukr = ukx * xr + uky * yr + ukz * zr;
429
430 // get the induced dipole field used for dipole torques
431
432 real tuir, tukr;
433
434 real tix3 = 2.0 * dscale * sr3 * ukx;
435 real tiy3 = 2.0 * dscale * sr3 * uky;
436 real tiz3 = 2.0 * dscale * sr3 * ukz;
437 real tkx3 = 2.0 * dscale * sr3 * uix;
438 real tky3 = 2.0 * dscale * sr3 * uiy;
439 real tkz3 = 2.0 * dscale * sr3 * uiz;
440 tuir = -2.0 * dscale * sr5 * ukr;
441 tukr = -2.0 * dscale * sr5 * uir;
442
443 uf0i += f * (tix3 + xr * tuir);
444 uf1i += f * (tiy3 + yr * tuir);
445 uf2i += f * (tiz3 + zr * tuir);
446 uf0k += f * (tkx3 + xr * tukr);
447 uf1k += f * (tky3 + yr * tukr);
448 uf2k += f * (tkz3 + zr * tukr);
449
450 // get induced dipole field gradient used for quadrupole torques
451
452 real tix5 = 4 * (dscale * sr5 * ukx);
453 real tiy5 = 4 * (dscale * sr5 * uky);
454 real tiz5 = 4 * (dscale * sr5 * ukz);
455 real tkx5 = 4 * (dscale * sr5 * uix);
456 real tky5 = 4 * (dscale * sr5 * uiy);
457 real tkz5 = 4 * (dscale * sr5 * uiz);
458 tuir = -2.0 * dscale * sr7 * ukr;
459 tukr = -2.0 * dscale * sr7 * uir;
460
461 duf0i += f * (xr * tix5 + xr * xr * tuir);
462 duf1i += f * (xr * tiy5 + yr * tix5 + 2 * xr * yr * tuir);
463 duf2i += f * (yr * tiy5 + yr * yr * tuir);
464 duf3i += f * (xr * tiz5 + zr * tix5 + 2 * xr * zr * tuir);
465 duf4i += f * (yr * tiz5 + zr * tiy5 + 2 * yr * zr * tuir);
466 duf5i += f * (zr * tiz5 + zr * zr * tuir);
467 duf0k += f * (-xr * tkx5 - xr * xr * tukr);
468 duf1k += f * (-xr * tky5 - yr * tkx5 - 2 * xr * yr * tukr);
469 duf2k += f * (-yr * tky5 - yr * yr * tukr);
470 duf3k += f * (-xr * tkz5 - zr * tkx5 - 2 * xr * zr * tukr);
471 duf4k += f * (-yr * tkz5 - zr * tky5 - 2 * yr * zr * tukr);
472 duf5k += f * (-zr * tkz5 - zr * zr * tukr);
473
474 // get the field gradient for direct polarization force
475
476 real term1, term2, term3, term4, term5, term6, term7;
477
478 term1 = bn[2] - ex3 * rr5;
479 term2 = bn[3] - ex5 * rr7;
480 term3 = -sr3 + term1 * xr * xr - xr * rc31;
481 term4 = rc31 - term1 * xr - sr5 * xr;
482 term5 = term2 * xr * xr - sr5 - xr * rc51;
483 term6 = (bn[4] - ex7 * rr9) * xr * xr - bn[3] - xr * rc71;
484 term7 = rc51 - 2 * bn[3] * xr + (ex5 + 1.5f * ex7) * rr7 * xr;
485 real tixx = ci * term3 + dix * term4 + dir * term5 + 2 * sr5 * qixx
486 + (qiy * yr + qiz * zr) * ex7 * rr7 + 2 * qix * term7 + qir * term6;
487 real tkxx = ck * term3 - dkx * term4 - dkr * term5 + 2 * sr5 * qkxx
488 + (qky * yr + qkz * zr) * ex7 * rr7 + 2 * qkx * term7 + qkr * term6;
489
490 term3 = -sr3 + term1 * yr * yr - yr * rc32;
491 term4 = rc32 - term1 * yr - sr5 * yr;
492 term5 = term2 * yr * yr - sr5 - yr * rc52;
493 term6 = (bn[4] - ex7 * rr9) * yr * yr - bn[3] - yr * rc72;
494 term7 = rc52 - 2 * bn[3] * yr + (ex5 + 1.5f * ex7) * rr7 * yr;
495 real tiyy = ci * term3 + diy * term4 + dir * term5 + 2 * sr5 * qiyy
496 + (qix * xr + qiz * zr) * ex7 * rr7 + 2 * qiy * term7 + qir * term6;
497 real tkyy = ck * term3 - dky * term4 - dkr * term5 + 2 * sr5 * qkyy
498 + (qkx * xr + qkz * zr) * ex7 * rr7 + 2 * qky * term7 + qkr * term6;
499
500 term3 = -sr3 + term1 * zr * zr - zr * rc33;
501 term4 = rc33 - term1 * zr - sr5 * zr;
502 term5 = term2 * zr * zr - sr5 - zr * rc53;
503 term6 = (bn[4] - ex7 * rr9) * zr * zr - bn[3] - zr * rc73;
504 term7 = rc53 - 2 * bn[3] * zr + (ex5 + 1.5f * ex7) * rr7 * zr;
505 real tizz = ci * term3 + diz * term4 + dir * term5 + 2 * sr5 * qizz
506 + (qix * xr + qiy * yr) * ex7 * rr7 + 2 * qiz * term7 + qir * term6;
507 real tkzz = ck * term3 - dkz * term4 - dkr * term5 + 2 * sr5 * qkzz
508 + (qkx * xr + qky * yr) * ex7 * rr7 + 2 * qkz * term7 + qkr * term6;
509
510 term3 = term1 * xr * yr - yr * rc31;
511 term4 = rc31 - term1 * xr;
512 term5 = term2 * xr * yr - yr * rc51;
513 term6 = (bn[4] - ex7 * rr9) * xr * yr - yr * rc71;
514 term7 = rc51 - term2 * xr;
515 real tixy = ci * term3 - sr5 * dix * yr + diy * term4 + dir * term5
516 + 2 * sr5 * qixy - 2 * sr7 * yr * qix + 2 * qiy * term7 + qir * term6;
517 real tkxy = ck * term3 + sr5 * dkx * yr - dky * term4 - dkr * term5
518 + 2 * sr5 * qkxy - 2 * sr7 * yr * qkx + 2 * qky * term7 + qkr * term6;
519
520 term3 = term1 * xr * zr - zr * rc31;
521 term5 = term2 * xr * zr - zr * rc51;
522 term6 = (bn[4] - ex7 * rr9) * xr * zr - zr * rc71;
523 real tixz = ci * term3 - sr5 * dix * zr + diz * term4 + dir * term5
524 + 2 * sr5 * qixz - 2 * sr7 * zr * qix + 2 * qiz * term7 + qir * term6;
525 real tkxz = ck * term3 + sr5 * dkx * zr - dkz * term4 - dkr * term5
526 + 2 * sr5 * qkxz - 2 * sr7 * zr * qkx + 2 * qkz * term7 + qkr * term6;
527
528 term3 = term1 * yr * zr - zr * rc32;
529 term4 = rc32 - term1 * yr;
530 term5 = term2 * yr * zr - zr * rc52;
531 term6 = (bn[4] - ex7 * rr9) * yr * zr - zr * rc72;
532 term7 = rc52 - term2 * yr;
533 real tiyz = ci * term3 - sr5 * diy * zr + diz * term4 + dir * term5
534 + 2 * sr5 * qiyz - 2 * sr7 * zr * qiy + 2 * qiz * term7 + qir * term6;
535 real tkyz = ck * term3 + sr5 * dky * zr - dkz * term4 - dkr * term5
536 + 2 * sr5 * qkyz - 2 * sr7 * zr * qky + 2 * qkz * term7 + qkr * term6;
537
538 // get the dEd/dR terms for Thole direct polarization force
539
540 real depx, depy, depz;
541 real frcx, frcy, frcz;
542
543 depx = tixx * ukx + tixy * uky + tixz * ukz - tkxx * uix - tkxy * uiy
544 - tkxz * uiz;
545 depy = tixy * ukx + tiyy * uky + tiyz * ukz - tkxy * uix - tkyy * uiy
546 - tkyz * uiz;
547 depz = tixz * ukx + tiyz * uky + tizz * ukz - tkxz * uix - tkyz * uiy
548 - tkzz * uiz;
549 if CONSTEXPR (eq<ETYP, EWALD>()) {
550 frcx = -2 * depx;
551 frcy = -2 * depy;
552 frcz = -2 * depz;
553 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
554 frcx = -2 * depx * dscale;
555 frcy = -2 * depy * dscale;
556 frcz = -2 * depz * dscale;
557 }
558
559 // get the dtau/dr terms used for mutual polarization force
560 damp_thole3g( //
561 r, rr2, xr, yr, zr, //
562 pdi, pti, pdk, ptk, //
563 ex3, ex5, ex7, //
564 rc31, rc32, rc33, //
565 rc51, rc52, rc53, //
566 rc71, rc72, rc73);
567 rc31 *= rr3;
568 rc32 *= rr3;
569 rc33 *= rr3;
570 rc51 *= rr5;
571 rc52 *= rr5;
572 rc53 *= rr5;
573 sr3 = bn[1] - ex3 * rr3;
574 sr5 = bn[2] - ex5 * rr5;
575
576 term1 = bn[2] - ex3 * rr5;
577 term2 = bn[3] - ex5 * rr7;
578 term3 = sr5 + term1;
579
580 term5 = -xr * term3 + rc31;
581 term6 = -sr5 + xr * xr * term2 - xr * rc51;
582 tixx = uix * term5 + uir * term6;
583 tkxx = ukx * term5 + ukr * term6;
584
585 term5 = -yr * term3 + rc32;
586 term6 = -sr5 + yr * yr * term2 - yr * rc52;
587 tiyy = uiy * term5 + uir * term6;
588 tkyy = uky * term5 + ukr * term6;
589
590 term5 = -zr * term3 + rc33;
591 term6 = -sr5 + zr * zr * term2 - zr * rc53;
592 tizz = uiz * term5 + uir * term6;
593 tkzz = ukz * term5 + ukr * term6;
594
595 term4 = -sr5 * yr;
596 term5 = -xr * term1 + rc31;
597 term6 = xr * yr * term2 - yr * rc51;
598 tixy = uix * term4 + uiy * term5 + uir * term6;
599 tkxy = ukx * term4 + uky * term5 + ukr * term6;
600
601 term4 = -sr5 * zr;
602 term6 = xr * zr * term2 - zr * rc51;
603 tixz = uix * term4 + uiz * term5 + uir * term6;
604 tkxz = ukx * term4 + ukz * term5 + ukr * term6;
605
606 term5 = -yr * term1 + rc32;
607 term6 = yr * zr * term2 - zr * rc52;
608 tiyz = uiy * term4 + uiz * term5 + uir * term6;
609 tkyz = uky * term4 + ukz * term5 + ukr * term6;
610
611 depx = tixx * ukx + tixy * uky + tixz * ukz + tkxx * uix + tkxy * uiy
612 + tkxz * uiz;
613 depy = tixy * ukx + tiyy * uky + tiyz * ukz + tkxy * uix + tkyy * uiy
614 + tkyz * uiz;
615 depz = tixz * ukx + tiyz * uky + tizz * ukz + tkxz * uix + tkyz * uiy
616 + tkzz * uiz;
617 if CONSTEXPR (eq<ETYP, EWALD>()) {
618 frcx -= depx;
619 frcy -= depy;
620 frcz -= depz;
621 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
622 frcx -= uscale * depx;
623 frcy -= uscale * depy;
624 frcz -= uscale * depz;
625 }
626
627 frcx *= f;
628 frcy *= f;
629 frcz *= f;
630 frcxi += frcx;
631 frcyi += frcy;
632 frczi += frcz;
633 frcxk -= frcx;
634 frcyk -= frcy;
635 frczk -= frcz;
636
637 if CONSTEXPR (do_v) {
638 vxx = -xr * frcx;
639 vxy = -0.5f * (yr * frcx + xr * frcy);
640 vxz = -0.5f * (zr * frcx + xr * frcz);
641 vyy = -yr * frcy;
642 vyz = -0.5f * (zr * frcy + yr * frcz);
643 vzz = -zr * frcz;
644 }
645 }
646}
647}
SEQ_CUDA
#define SEQ_CUDA
Definition: acc/seqdef.h:12
MAYBE_UNUSED
#define MAYBE_UNUSED
Definition: macro.h:75
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::pair_polar_aplus_v2
__device__ void pair_polar_aplus_v2(real r2, real xr, real yr, real zr, real dscale, real uscale, real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real uix, real uiy, real uiz, real pdi, real pti, real ddi, real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real ukx, real uky, real ukz, real pdk, real ptk, real ddk, real f, real aewald, real &__restrict__ frcxi, real &__restrict__ frcyi, real &__restrict__ frczi, real &__restrict__ frcxk, real &__restrict__ frcyk, real &__restrict__ frczk, real &__restrict__ uf0i, real &__restrict__ uf1i, real &__restrict__ uf2i, real &__restrict__ uf0k, real &__restrict__ uf1k, real &__restrict__ uf2k, real &__restrict__ duf0i, real &__restrict__ duf1i, real &__restrict__ duf2i, real &__restrict__ duf3i, real &__restrict__ duf4i, real &__restrict__ duf5i, real &__restrict__ duf0k, real &__restrict__ duf1k, real &__restrict__ duf2k, real &__restrict__ duf3k, real &__restrict__ duf4k, real &__restrict__ duf5k, real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vxy, real &__restrict__ vxz, real &__restrict__ vyy, real &__restrict__ vyz, real &__restrict__ vzz, real &__restrict__ poti, real &__restrict__ potk)
Definition: pairpolaraplus.h:312
tinker::depy
grad_prec * depy
tinker::depx
grad_prec * depx
tinker::damp_thole3g
__device__ void damp_thole3g(real r, real rr2, real xr, real yr, real zr, real pdi, real pti, real pdk, real ptk, real &__restrict__ scale31, real &__restrict__ scale51, real &__restrict__ scale71, real &__restrict__ rc31, real &__restrict__ rc32, real &__restrict__ rc33, real &__restrict__ rc51, real &__restrict__ rc52, real &__restrict__ rc53, real &__restrict__ rc71, real &__restrict__ rc72, real &__restrict__ rc73)
Definition: damp.h:79
tinker::damp_aplus3g
__device__ void damp_aplus3g(real r, real rr2, real xr, real yr, real zr, real pdi, real ddi, real pdk, real ddk, real &__restrict__ scale31, real &__restrict__ scale51, real &__restrict__ scale71, real &__restrict__ rc31, real &__restrict__ rc32, real &__restrict__ rc33, real &__restrict__ rc51, real &__restrict__ rc52, real &__restrict__ rc53, real &__restrict__ rc71, real &__restrict__ rc72, real &__restrict__ rc73)
Definition: dampaplus.h:146
tinker::depz
grad_prec * depz
tinker::damp_aplus3
__device__ void damp_aplus3(real r, real pdi, real ddi, real pdk, real ddk, real &__restrict__ scale3, real &__restrict__ scale5, real &__restrict__ scale7)
Definition: dampaplus.h:131
tinker::pair_polar_aplus
__device__ void pair_polar_aplus(real r2, real xr, real yr, real zr, real dscale, real uscale, real ci, real dix, real diy, real diz, real pdi, real pti, real ddi, 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 pdk, real ptk, real ddk, 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: pairpolaraplus.h:10
tinker::PairPolarGrad
Definition: pair_polar.h:8
Generated by doxygen 1.9.5