Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pair_mpole.h
1#pragma once
2#include "ff/elec.h"
3#include "seq/add.h"
4#include "seq/damp.h"
5
6namespace tinker {
13struct PairMPoleGrad
14{
18 real frcx, frcy, frcz;
21 real ttmi[3];
23 real ttmk[3];
24};
25
26SEQ_ROUTINE
27inline void zero(PairMPoleGrad& pgrad)
28{
29 pgrad.frcx = 0;
30 pgrad.frcy = 0;
31 pgrad.frcz = 0;
32 pgrad.ttmi[0] = 0;
33 pgrad.ttmi[1] = 0;
34 pgrad.ttmi[2] = 0;
35 pgrad.ttmk[0] = 0;
36 pgrad.ttmk[1] = 0;
37 pgrad.ttmk[2] = 0;
38}
39
62#pragma acc routine seq
63template <bool do_e, bool do_g, class ETYP>
64SEQ_CUDA
65void pair_mpole( //
66 real r2, real xr, real yr, real zr, real mscale, //
67 real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz,
68 real qiyy, real qiyz,
69 real qizz, //
70 real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz,
71 real qkyy, real qkyz,
72 real qkzz, //
73 real f, real aewald, real& restrict e, PairMPoleGrad& restrict pgrad)
74{
75 real r = REAL_SQRT(r2);
76 real invr1 = REAL_RECIP(r);
77 real rr2 = invr1 * invr1;
78
79 real bn[6];
80 real& rr1 = bn[0];
81 real& rr3 = bn[1];
82 real& rr5 = bn[2];
83 real& rr7 = bn[3];
84 real& rr9 = bn[4];
85 real& rr11 = bn[5];
86
87 if CONSTEXPR (eq<ETYP, EWALD>()) {
88 if CONSTEXPR (!do_g)
89 damp_ewald<5>(bn, r, invr1, rr2, aewald);
90 else
91 damp_ewald<6>(bn, r, invr1, rr2, aewald);
92 bn[0] *= f;
93 bn[1] *= f;
94 bn[2] *= f;
95 bn[3] *= f;
96 bn[4] *= f;
97 if CONSTEXPR (do_g) bn[5] *= f;
98 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
99 rr1 = mscale * f * invr1;
100 rr3 = rr1 * rr2;
101 rr5 = 3 * rr3 * rr2;
102 rr7 = 5 * rr5 * rr2;
103 rr9 = 7 * rr7 * rr2;
104 if CONSTEXPR (do_g) rr11 = 9 * rr9 * rr2;
105 }
106
107 real dir = dix * xr + diy * yr + diz * zr;
108 real qix = qixx * xr + qixy * yr + qixz * zr;
109 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
110 real qiz = qixz * xr + qiyz * yr + qizz * zr;
111 real qir = qix * xr + qiy * yr + qiz * zr;
112 real dkr = dkx * xr + dky * yr + dkz * zr;
113 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
114 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
115 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
116 real qkr = qkx * xr + qky * yr + qkz * zr;
117 real dik = dix * dkx + diy * dky + diz * dkz;
118 real qik = qix * qkx + qiy * qky + qiz * qkz;
119 real diqk = dix * qkx + diy * qky + diz * qkz;
120 real dkqi = dkx * qix + dky * qiy + dkz * qiz;
121 real qiqk = 2 * (qixy * qkxy + qixz * qkxz + qiyz * qkyz) + qixx * qkxx
122 + qiyy * qkyy + qizz * qkzz;
123
124 real term1 = ci * ck;
125 real term2 = ck * dir - ci * dkr + dik;
126 real term3 = ci * qkr + ck * qir - dir * dkr + 2 * (dkqi - diqk + qiqk);
127 real term4 = dir * qkr - dkr * qir - 4 * qik;
128 real term5 = qir * qkr;
129
130 if CONSTEXPR (do_e) {
131 e = term1 * rr1 + term2 * rr3 + term3 * rr5 + term4 * rr7 + term5 * rr9;
132 } // end if (do_e)
133
134 if CONSTEXPR (do_g) {
135
136 // gradient
137
138 real qixk = qixx * qkx + qixy * qky + qixz * qkz;
139 real qiyk = qixy * qkx + qiyy * qky + qiyz * qkz;
140 real qizk = qixz * qkx + qiyz * qky + qizz * qkz;
141 real qkxi = qkxx * qix + qkxy * qiy + qkxz * qiz;
142 real qkyi = qkxy * qix + qkyy * qiy + qkyz * qiz;
143 real qkzi = qkxz * qix + qkyz * qiy + qkzz * qiz;
144
145 real diqkx = dix * qkxx + diy * qkxy + diz * qkxz;
146 real diqky = dix * qkxy + diy * qkyy + diz * qkyz;
147 real diqkz = dix * qkxz + diy * qkyz + diz * qkzz;
148 real dkqix = dkx * qixx + dky * qixy + dkz * qixz;
149 real dkqiy = dkx * qixy + dky * qiyy + dkz * qiyz;
150 real dkqiz = dkx * qixz + dky * qiyz + dkz * qizz;
151
152 real de = term1 * rr3 + term2 * rr5 + term3 * rr7 + term4 * rr9
153 + term5 * rr11;
154
155 term1 = -ck * rr3 + dkr * rr5 - qkr * rr7;
156 term2 = ci * rr3 + dir * rr5 + qir * rr7;
157 term3 = 2 * rr5;
158 term4 = 2 * (-ck * rr5 + dkr * rr7 - qkr * rr9);
159 term5 = 2 * (-ci * rr5 - dir * rr7 - qir * rr9);
160 real term6 = 4 * rr7;
161
162 pgrad.frcx = de * xr + term1 * dix + term2 * dkx + term3 * (diqkx - dkqix)
163 + term4 * qix + term5 * qkx + term6 * (qixk + qkxi);
164 pgrad.frcy = de * yr + term1 * diy + term2 * dky + term3 * (diqky - dkqiy)
165 + term4 * qiy + term5 * qky + term6 * (qiyk + qkyi);
166 pgrad.frcz = de * zr + term1 * diz + term2 * dkz + term3 * (diqkz - dkqiz)
167 + term4 * qiz + term5 * qkz + term6 * (qizk + qkzi);
168
169 // torque
170
171 real dirx = diy * zr - diz * yr;
172 real diry = diz * xr - dix * zr;
173 real dirz = dix * yr - diy * xr;
174 real dkrx = dky * zr - dkz * yr;
175 real dkry = dkz * xr - dkx * zr;
176 real dkrz = dkx * yr - dky * xr;
177 real dikx = diy * dkz - diz * dky;
178 real diky = diz * dkx - dix * dkz;
179 real dikz = dix * dky - diy * dkx;
180
181 real qirx = qiz * yr - qiy * zr;
182 real qiry = qix * zr - qiz * xr;
183 real qirz = qiy * xr - qix * yr;
184 real qkrx = qkz * yr - qky * zr;
185 real qkry = qkx * zr - qkz * xr;
186 real qkrz = qky * xr - qkx * yr;
187 real qikx = qky * qiz - qkz * qiy;
188 real qiky = qkz * qix - qkx * qiz;
189 real qikz = qkx * qiy - qky * qix;
190
191 real qikrx = qizk * yr - qiyk * zr;
192 real qikry = qixk * zr - qizk * xr;
193 real qikrz = qiyk * xr - qixk * yr;
194 real qkirx = qkzi * yr - qkyi * zr;
195 real qkiry = qkxi * zr - qkzi * xr;
196 real qkirz = qkyi * xr - qkxi * yr;
197
198 real diqkrx = diqkz * yr - diqky * zr;
199 real diqkry = diqkx * zr - diqkz * xr;
200 real diqkrz = diqky * xr - diqkx * yr;
201 real dkqirx = dkqiz * yr - dkqiy * zr;
202 real dkqiry = dkqix * zr - dkqiz * xr;
203 real dkqirz = dkqiy * xr - dkqix * yr;
204
205 real dqikx = diy * qkz - diz * qky + dky * qiz - dkz * qiy
206 - 2
207 * (qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy
208 - qiyz * qkyy - qizz * qkyz);
209 real dqiky = diz * qkx - dix * qkz + dkz * qix - dkx * qiz
210 - 2
211 * (qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz
212 - qixy * qkyz - qixz * qkzz);
213 real dqikz = dix * qky - diy * qkx + dkx * qiy - dky * qix
214 - 2
215 * (qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx
216 - qiyy * qkxy - qiyz * qkxz);
217
218 pgrad.ttmi[0] = -rr3 * dikx + term1 * dirx + term3 * (dqikx + dkqirx)
219 - term4 * qirx - term6 * (qikrx + qikx);
220 pgrad.ttmi[1] = -rr3 * diky + term1 * diry + term3 * (dqiky + dkqiry)
221 - term4 * qiry - term6 * (qikry + qiky);
222 pgrad.ttmi[2] = -rr3 * dikz + term1 * dirz + term3 * (dqikz + dkqirz)
223 - term4 * qirz - term6 * (qikrz + qikz);
224 pgrad.ttmk[0] = rr3 * dikx + term2 * dkrx - term3 * (dqikx + diqkrx)
225 - term5 * qkrx - term6 * (qkirx - qikx);
226 pgrad.ttmk[1] = rr3 * diky + term2 * dkry - term3 * (dqiky + diqkry)
227 - term5 * qkry - term6 * (qkiry - qiky);
228 pgrad.ttmk[2] = rr3 * dikz + term2 * dkrz - term3 * (dqikz + diqkrz)
229 - term5 * qkrz - term6 * (qkirz - qikz);
230 } // end if (do_g)
231}
232
233#pragma acc routine seq
234template <class Ver, class ETYP>
235SEQ_CUDA
236void pair_mpole_v2(real r2, real xr, real yr, real zr, real mscale, //
237 real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz,
238 real qiyy, real qiyz,
239 real qizz, //
240 real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz,
241 real qkyy, real qkyz,
242 real qkzz, //
243 real f, real aewald, //
244 real& restrict frcxi, real& restrict frcyi, real& restrict frczi,
245 real& restrict frcxk, real& restrict frcyk, real& restrict frczk,
246 real& restrict trqxi, real& restrict trqyi, real& restrict trqzi,
247 real& restrict trqxk, real& restrict trqyk, real& restrict trqzk, //
248 real& restrict e, //
249 real& restrict vxx, real& restrict vxy, real& restrict vxz,
250 real& restrict vyy, real& restrict vyz, real& restrict vzz)
251{
252 constexpr bool do_e = Ver::e;
253 constexpr bool do_g = Ver::g;
254 constexpr bool do_v = Ver::v;
255
256 real r = REAL_SQRT(r2);
257 real invr1 = REAL_RECIP(r);
258 real rr2 = invr1 * invr1;
259
260 real bn[6];
261 real& rr1 = bn[0];
262 real& rr3 = bn[1];
263 real& rr5 = bn[2];
264 real& rr7 = bn[3];
265 real& rr9 = bn[4];
266 real& rr11 = bn[5];
267
268 if CONSTEXPR (eq<ETYP, EWALD>()) {
269 if CONSTEXPR (!do_g)
270 damp_ewald<5>(bn, r, invr1, rr2, aewald);
271 else
272 damp_ewald<6>(bn, r, invr1, rr2, aewald);
273 bn[0] *= f;
274 bn[1] *= f;
275 bn[2] *= f;
276 bn[3] *= f;
277 bn[4] *= f;
278 if CONSTEXPR (do_g) bn[5] *= f;
279 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
280 rr1 = mscale * f * invr1;
281 rr3 = rr1 * rr2;
282 rr5 = 3 * rr3 * rr2;
283 rr7 = 5 * rr5 * rr2;
284 rr9 = 7 * rr7 * rr2;
285 if CONSTEXPR (do_g) rr11 = 9 * rr9 * rr2;
286 }
287
288 real dir = dix * xr + diy * yr + diz * zr;
289 real qix = qixx * xr + qixy * yr + qixz * zr;
290 real qiy = qixy * xr + qiyy * yr + qiyz * zr;
291 real qiz = qixz * xr + qiyz * yr + qizz * zr;
292 real qir = qix * xr + qiy * yr + qiz * zr;
293 real dkr = dkx * xr + dky * yr + dkz * zr;
294 real qkx = qkxx * xr + qkxy * yr + qkxz * zr;
295 real qky = qkxy * xr + qkyy * yr + qkyz * zr;
296 real qkz = qkxz * xr + qkyz * yr + qkzz * zr;
297 real qkr = qkx * xr + qky * yr + qkz * zr;
298 real dik = dix * dkx + diy * dky + diz * dkz;
299 real qik = qix * qkx + qiy * qky + qiz * qkz;
300 real diqk = dix * qkx + diy * qky + diz * qkz;
301 real dkqi = dkx * qix + dky * qiy + dkz * qiz;
302 real qiqk = 2 * (qixy * qkxy + qixz * qkxz + qiyz * qkyz) + qixx * qkxx
303 + qiyy * qkyy + qizz * qkzz;
304
305 real term1 = ci * ck;
306 real term2 = ck * dir - ci * dkr + dik;
307 real term3 = ci * qkr + ck * qir - dir * dkr + 2 * (dkqi - diqk + qiqk);
308 real term4 = dir * qkr - dkr * qir - 4 * qik;
309 real term5 = qir * qkr;
310
311 if CONSTEXPR (do_e) {
312 e = term1 * rr1 + term2 * rr3 + term3 * rr5 + term4 * rr7 + term5 * rr9;
313 } // end if (do_e)
314
315 real frcx, frcy, frcz;
316 if CONSTEXPR (do_g) {
317 // gradient
318 real qixk = qixx * qkx + qixy * qky + qixz * qkz;
319 real qiyk = qixy * qkx + qiyy * qky + qiyz * qkz;
320 real qizk = qixz * qkx + qiyz * qky + qizz * qkz;
321 real qkxi = qkxx * qix + qkxy * qiy + qkxz * qiz;
322 real qkyi = qkxy * qix + qkyy * qiy + qkyz * qiz;
323 real qkzi = qkxz * qix + qkyz * qiy + qkzz * qiz;
324
325 real diqkx = dix * qkxx + diy * qkxy + diz * qkxz;
326 real diqky = dix * qkxy + diy * qkyy + diz * qkyz;
327 real diqkz = dix * qkxz + diy * qkyz + diz * qkzz;
328 real dkqix = dkx * qixx + dky * qixy + dkz * qixz;
329 real dkqiy = dkx * qixy + dky * qiyy + dkz * qiyz;
330 real dkqiz = dkx * qixz + dky * qiyz + dkz * qizz;
331
332 real de = term1 * rr3 + term2 * rr5 + term3 * rr7 + term4 * rr9
333 + term5 * rr11;
334
335 term1 = -ck * rr3 + dkr * rr5 - qkr * rr7;
336 term2 = ci * rr3 + dir * rr5 + qir * rr7;
337 term3 = 2 * rr5;
338 term4 = 2 * (-ck * rr5 + dkr * rr7 - qkr * rr9);
339 term5 = 2 * (-ci * rr5 - dir * rr7 - qir * rr9);
340 real term6 = 4 * rr7;
341
342 frcx = de * xr + term1 * dix + term2 * dkx + term3 * (diqkx - dkqix)
343 + term4 * qix + term5 * qkx + term6 * (qixk + qkxi);
344 frcy = de * yr + term1 * diy + term2 * dky + term3 * (diqky - dkqiy)
345 + term4 * qiy + term5 * qky + term6 * (qiyk + qkyi);
346 frcz = de * zr + term1 * diz + term2 * dkz + term3 * (diqkz - dkqiz)
347 + term4 * qiz + term5 * qkz + term6 * (qizk + qkzi);
348 frcxi += frcx;
349 frcyi += frcy;
350 frczi += frcz;
351 frcxk -= frcx;
352 frcyk -= frcy;
353 frczk -= frcz;
354
355 // torque
356 real dirx = diy * zr - diz * yr;
357 real diry = diz * xr - dix * zr;
358 real dirz = dix * yr - diy * xr;
359 real dkrx = dky * zr - dkz * yr;
360 real dkry = dkz * xr - dkx * zr;
361 real dkrz = dkx * yr - dky * xr;
362 real dikx = diy * dkz - diz * dky;
363 real diky = diz * dkx - dix * dkz;
364 real dikz = dix * dky - diy * dkx;
365
366 real qirx = qiz * yr - qiy * zr;
367 real qiry = qix * zr - qiz * xr;
368 real qirz = qiy * xr - qix * yr;
369 real qkrx = qkz * yr - qky * zr;
370 real qkry = qkx * zr - qkz * xr;
371 real qkrz = qky * xr - qkx * yr;
372 real qikx = qky * qiz - qkz * qiy;
373 real qiky = qkz * qix - qkx * qiz;
374 real qikz = qkx * qiy - qky * qix;
375
376 real qikrx = qizk * yr - qiyk * zr;
377 real qikry = qixk * zr - qizk * xr;
378 real qikrz = qiyk * xr - qixk * yr;
379 real qkirx = qkzi * yr - qkyi * zr;
380 real qkiry = qkxi * zr - qkzi * xr;
381 real qkirz = qkyi * xr - qkxi * yr;
382
383 real diqkrx = diqkz * yr - diqky * zr;
384 real diqkry = diqkx * zr - diqkz * xr;
385 real diqkrz = diqky * xr - diqkx * yr;
386 real dkqirx = dkqiz * yr - dkqiy * zr;
387 real dkqiry = dkqix * zr - dkqiz * xr;
388 real dkqirz = dkqiy * xr - dkqix * yr;
389
390 real dqikx = diy * qkz - diz * qky + dky * qiz - dkz * qiy
391 - 2
392 * (qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy
393 - qiyz * qkyy - qizz * qkyz);
394 real dqiky = diz * qkx - dix * qkz + dkz * qix - dkx * qiz
395 - 2
396 * (qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz
397 - qixy * qkyz - qixz * qkzz);
398 real dqikz = dix * qky - diy * qkx + dkx * qiy - dky * qix
399 - 2
400 * (qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx
401 - qiyy * qkxy - qiyz * qkxz);
402
403 trqxi += -rr3 * dikx + term1 * dirx + term3 * (dqikx + dkqirx)
404 - term4 * qirx - term6 * (qikrx + qikx);
405 trqyi += -rr3 * diky + term1 * diry + term3 * (dqiky + dkqiry)
406 - term4 * qiry - term6 * (qikry + qiky);
407 trqzi += -rr3 * dikz + term1 * dirz + term3 * (dqikz + dkqirz)
408 - term4 * qirz - term6 * (qikrz + qikz);
409 trqxk += rr3 * dikx + term2 * dkrx - term3 * (dqikx + diqkrx)
410 - term5 * qkrx - term6 * (qkirx - qikx);
411 trqyk += rr3 * diky + term2 * dkry - term3 * (dqiky + diqkry)
412 - term5 * qkry - term6 * (qkiry - qiky);
413 trqzk += rr3 * dikz + term2 * dkrz - term3 * (dqikz + diqkrz)
414 - term5 * qkrz - term6 * (qkirz - qikz);
415 } // end if (do_g)
416
417 if CONSTEXPR (do_v) {
418 vxx = -xr * frcx;
419 vxy = -0.5f * (yr * frcx + xr * frcy);
420 vxz = -0.5f * (zr * frcx + xr * frcz);
421 vyy = -yr * frcy;
422 vyz = -0.5f * (zr * frcy + yr * frcz);
423 vzz = -zr * frcz;
424 } // end if (do_v)
425}
426}
SEQ_ROUTINE
#define SEQ_ROUTINE
Definition: acc/seqdef.h:7
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::PairMPoleGrad::frcy
real frcy
Incomplete x, y, or z gradients of atom i.
Definition: pair_mpole.h:18
tinker::PairMPoleGrad::frcx
real frcx
Incomplete x, y, or z gradients of atom i.
Definition: pair_mpole.h:18
tinker::PairMPoleGrad::ttmi
real ttmi[3]
x, y, and z torques on atom i.
Definition: pair_mpole.h:21
tinker::PairMPoleGrad::frcz
real frcz
Incomplete x, y, or z gradients of atom i.
Definition: pair_mpole.h:18
tinker::PairMPoleGrad::ttmk
real ttmk[3]
x, y, and z torques on atom k.
Definition: pair_mpole.h:23
tinker::pair_mpole
__device__ void pair_mpole(real r2, real xr, real yr, real zr, real mscale, real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real f, real aewald, real &__restrict__ e, PairMPoleGrad &__restrict__ pgrad)
OpenACC pairwise multipole electrostatic energy.
Definition: pair_mpole.h:65
tinker::PairMPoleGrad
Components of the pairwise multipole electrostatic gradient between atoms i and k,...
Definition: pair_mpole.h:14
tinker::real
float real
Definition: precision.h:80
tinker
Definition: testrt.h:9
tinker::zero
__device__ void zero(PairMPoleGrad &pgrad)
Definition: pair_mpole.h:27
tinker::pair_mpole_v2
__device__ void pair_mpole_v2(real r2, real xr, real yr, real zr, real mscale, real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real f, real aewald, real &__restrict__ frcxi, real &__restrict__ frcyi, real &__restrict__ frczi, real &__restrict__ frcxk, real &__restrict__ frcyk, real &__restrict__ frczk, real &__restrict__ trqxi, real &__restrict__ trqyi, real &__restrict__ trqzi, real &__restrict__ trqxk, real &__restrict__ trqyk, real &__restrict__ trqzk, real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vxy, real &__restrict__ vxz, real &__restrict__ vyy, real &__restrict__ vyz, real &__restrict__ vzz)
Definition: pair_mpole.h:236
Generated by doxygen 1.9.5