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