Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pairmpoleaplus.h
1#pragma once
2#include "seq/dampaplus.h"
3#include "seq/pair_mpole.h"
4
5namespace tinker {
6#pragma acc routine seq
7template <bool do_e, bool do_g, class ETYP, int CFLX>
8SEQ_CUDA
9void pair_mpole_chgpen_aplus( //
10 real r2, real xr, real yr, real zr, real mscale, //
11 real ci, real dix, real diy, real diz, real corei, real vali, real alphai,
12 real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, //
13 real ck, real dkx, real dky, real dkz, real corek, real valk, real alphak,
14 real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, //
15 real f, real aewald, real& restrict e, real& restrict poti,
16 real& restrict potk, PairMPoleGrad& restrict pgrad)
17{
18 real r = REAL_SQRT(r2);
19 real invr1 = REAL_RECIP(r);
20 real rr2 = invr1 * invr1;
21
22 real dmpi[5];
23 real dmpk[5];
24 real dmpik[6];
25 real bn[6];
26
27 real rr1 = f * invr1;
28 real rr3 = rr1 * rr2;
29 real rr5 = 3 * rr3 * rr2;
30 real rr7 = 5 * rr5 * rr2;
31 real rr9 = 7 * rr7 * rr2;
32 real rr11;
33
34 if CONSTEXPR (do_g) rr11 = 9 * rr9 * 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 real dik = dix * dkx + diy * dky + diz * dkz;
47 real qik = qix * qkx + qiy * qky + qiz * qkz;
48 real diqk = dix * qkx + diy * qky + diz * qkz;
49 real dkqi = dkx * qix + dky * qiy + dkz * qiz;
50 real qiqk = 2 * (qixy * qkxy + qixz * qkxz + qiyz * qkyz) + qixx * qkxx
51 + qiyy * qkyy + qizz * qkzz;
52
53 // chgpen_aplus terms
54 real term1 = corei * corek;
55 real term1i = corek * vali;
56 real term2i = corek * dir;
57 real term3i = corek * qir;
58 real term1k = corei * valk;
59 real term2k = -corei * dkr;
60 real term3k = corei * qkr;
61 real term1ik = vali * valk;
62 real term2ik = valk * dir - vali * dkr + dik;
63 real term3ik = vali * qkr + valk * qir - dir * dkr
64 + 2 * (dkqi - diqk + qiqk);
65 real term4ik = dir * qkr - dkr * qir - 4 * qik;
66 real term5ik = qir * qkr;
67
68 real rr1i, rr3i, rr5i, rr7i, rr1k, rr3k, rr5k, rr7k, rr1ik, rr3ik, rr5ik,
69 rr7ik, rr9ik, rr11ik;
70
71 // Compute damping factors
72 if CONSTEXPR (do_g) {
73 damp_gordon2<11>(dmpik, dmpi, dmpk, r, alphai, alphak);
74 } else {
75 damp_gordon2<9>(dmpik, dmpi, dmpk, r, alphai, alphak);
76 }
77 //
78
79 if CONSTEXPR (eq<ETYP, EWALD>()) {
80 if CONSTEXPR (do_g) {
81 damp_ewald<6>(bn, r, invr1, rr2, aewald);
82 } else {
83 damp_ewald<5>(bn, r, invr1, rr2, aewald);
84 }
85
86 bn[0] *= f;
87 bn[1] *= f;
88 bn[2] *= f;
89 bn[3] *= f;
90 bn[4] *= f;
91 if CONSTEXPR (do_g) bn[5] *= f;
92 } else if CONSTEXPR (eq<ETYP, NON_EWALD>()) {
93 bn[0] = rr1;
94 bn[1] = rr3;
95 bn[2] = rr5;
96 bn[3] = rr7;
97 bn[4] = rr9;
98 if CONSTEXPR (do_g) bn[5] = rr11;
99 } // endif NON_EWALD
100
101 real m = 1 - mscale;
102 rr1i = bn[0] - (m + mscale * dmpi[0]) * rr1;
103 rr3i = bn[1] - (m + mscale * dmpi[1]) * rr3;
104 rr5i = bn[2] - (m + mscale * dmpi[2]) * rr5;
105 rr7i = bn[3] - (m + mscale * dmpi[3]) * rr7;
106 rr1k = bn[0] - (m + mscale * dmpk[0]) * rr1;
107 rr3k = bn[1] - (m + mscale * dmpk[1]) * rr3;
108 rr5k = bn[2] - (m + mscale * dmpk[2]) * rr5;
109 rr7k = bn[3] - (m + mscale * dmpk[3]) * rr7;
110 rr1ik = bn[0] - (m + mscale * dmpik[0]) * rr1;
111 rr3ik = bn[1] - (m + mscale * dmpik[1]) * rr3;
112 rr5ik = bn[2] - (m + mscale * dmpik[2]) * rr5;
113 rr7ik = bn[3] - (m + mscale * dmpik[3]) * rr7;
114 rr9ik = bn[4] - (m + mscale * dmpik[4]) * rr9;
115 rr1 = bn[0] - m * rr1;
116 rr3 = bn[1] - m * rr3;
117
118 if CONSTEXPR (do_g) rr11ik = bn[5] - (m + mscale * dmpik[5]) * rr11;
119
120 if CONSTEXPR (do_e)
121 e = term1 * rr1 + term4ik * rr7ik + term5ik * rr9ik + term1i * rr1i
122 + term1k * rr1k + term1ik * rr1ik + term2i * rr3i + term2k * rr3k
123 + term2ik * rr3ik + term3i * rr5i + term3k * rr5k + term3ik * rr5ik;
124 // end if (do_e)
125
126 if CONSTEXPR (do_g) {
127 // gradient
128 real qixk = qixx * qkx + qixy * qky + qixz * qkz;
129 real qiyk = qixy * qkx + qiyy * qky + qiyz * qkz;
130 real qizk = qixz * qkx + qiyz * qky + qizz * qkz;
131 real qkxi = qkxx * qix + qkxy * qiy + qkxz * qiz;
132 real qkyi = qkxy * qix + qkyy * qiy + qkyz * qiz;
133 real qkzi = qkxz * qix + qkyz * qiy + qkzz * qiz;
134
135 real diqkx = dix * qkxx + diy * qkxy + diz * qkxz;
136 real diqky = dix * qkxy + diy * qkyy + diz * qkyz;
137 real diqkz = dix * qkxz + diy * qkyz + diz * qkzz;
138 real dkqix = dkx * qixx + dky * qixy + dkz * qixz;
139 real dkqiy = dkx * qixy + dky * qiyy + dkz * qiyz;
140 real dkqiz = dkx * qixz + dky * qiyz + dkz * qizz;
141
142 real de = term1 * rr3 + term4ik * rr9ik + term5ik * rr11ik + term1i * rr3i
143 + term1k * rr3k + term1ik * rr3ik + term2i * rr5i + term2k * rr5k
144 + term2ik * rr5ik + term3i * rr7i + term3k * rr7k + term3ik * rr7ik;
145
146 term1 = -corek * rr3i - valk * rr3ik + dkr * rr5ik - qkr * rr7ik;
147 real term2 = corei * rr3k + vali * rr3ik + dir * rr5ik + qir * rr7ik;
148 real term3 = 2 * rr5ik;
149 real term4 = -2
150 * (corek * rr5i + valk * rr5ik - dkr * rr7ik + qkr * rr9ik);
151 real term5 = -2
152 * (corei * rr5k + vali * rr5ik + dir * rr7ik + qir * rr9ik);
153 real term6 = 4 * rr7ik;
154
155 if CONSTEXPR (CFLX) {
156 real t1i = corek * rr1i + valk * rr1ik;
157 real t1k = corei * rr1k + vali * rr1ik;
158 real t2i = -dkr * rr3ik;
159 real t2k = dir * rr3ik;
160 real t3i = qkr * rr5ik;
161 real t3k = qir * rr5ik;
162 poti = t1i + t2i + t3i;
163 potk = t1k + t2k + t3k;
164 }
165
166 pgrad.frcx = de * xr + term1 * dix + term2 * dkx + term3 * (diqkx - dkqix)
167 + term4 * qix + term5 * qkx + term6 * (qixk + qkxi);
168
169 pgrad.frcy = de * yr + term1 * diy + term2 * dky + term3 * (diqky - dkqiy)
170 + term4 * qiy + term5 * qky + term6 * (qiyk + qkyi);
171 pgrad.frcz = de * zr + term1 * diz + term2 * dkz + term3 * (diqkz - dkqiz)
172 + term4 * qiz + term5 * qkz + term6 * (qizk + qkzi);
173
174 // torque
175 real dirx = diy * zr - diz * yr;
176 real diry = diz * xr - dix * zr;
177 real dirz = dix * yr - diy * xr;
178 real dkrx = dky * zr - dkz * yr;
179 real dkry = dkz * xr - dkx * zr;
180 real dkrz = dkx * yr - dky * xr;
181 real dikx = diy * dkz - diz * dky;
182 real diky = diz * dkx - dix * dkz;
183 real dikz = dix * dky - diy * dkx;
184
185 real qirx = qiz * yr - qiy * zr;
186 real qiry = qix * zr - qiz * xr;
187 real qirz = qiy * xr - qix * yr;
188 real qkrx = qkz * yr - qky * zr;
189 real qkry = qkx * zr - qkz * xr;
190 real qkrz = qky * xr - qkx * yr;
191 real qikx = qky * qiz - qkz * qiy;
192 real qiky = qkz * qix - qkx * qiz;
193 real qikz = qkx * qiy - qky * qix;
194
195 real qikrx = qizk * yr - qiyk * zr;
196 real qikry = qixk * zr - qizk * xr;
197 real qikrz = qiyk * xr - qixk * yr;
198 real qkirx = qkzi * yr - qkyi * zr;
199 real qkiry = qkxi * zr - qkzi * xr;
200 real qkirz = qkyi * xr - qkxi * yr;
201
202 real diqkrx = diqkz * yr - diqky * zr;
203 real diqkry = diqkx * zr - diqkz * xr;
204 real diqkrz = diqky * xr - diqkx * yr;
205 real dkqirx = dkqiz * yr - dkqiy * zr;
206 real dkqiry = dkqix * zr - dkqiz * xr;
207 real dkqirz = dkqiy * xr - dkqix * yr;
208
209 real dqikx = diy * qkz - diz * qky + dky * qiz - dkz * qiy
210 - 2
211 * (qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy
212 - qiyz * qkyy - qizz * qkyz);
213 real dqiky = diz * qkx - dix * qkz + dkz * qix - dkx * qiz
214 - 2
215 * (qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz
216 - qixy * qkyz - qixz * qkzz);
217 real dqikz = dix * qky - diy * qkx + dkx * qiy - dky * qix
218 - 2
219 * (qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx
220 - qiyy * qkxy - qiyz * qkxz);
221
222 pgrad.ttmi[0] = -rr3ik * dikx + term1 * dirx + term3 * (dqikx + dkqirx)
223 - term4 * qirx - term6 * (qikrx + qikx);
224 pgrad.ttmi[1] = -rr3ik * diky + term1 * diry + term3 * (dqiky + dkqiry)
225 - term4 * qiry - term6 * (qikry + qiky);
226 pgrad.ttmi[2] = -rr3ik * dikz + term1 * dirz + term3 * (dqikz + dkqirz)
227 - term4 * qirz - term6 * (qikrz + qikz);
228 pgrad.ttmk[0] = rr3ik * dikx + term2 * dkrx - term3 * (dqikx + diqkrx)
229 - term5 * qkrx - term6 * (qkirx - qikx);
230 pgrad.ttmk[1] = rr3ik * diky + term2 * dkry - term3 * (dqiky + diqkry)
231 - term5 * qkry - term6 * (qkiry - qiky);
232 pgrad.ttmk[2] = rr3ik * dikz + term2 * dkrz - term3 * (dqikz + diqkrz)
233 - term5 * qkrz - term6 * (qkirz - qikz);
234 } // end if (do_g)
235}
236}
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
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::pair_mpole_chgpen_aplus
__device__ void pair_mpole_chgpen_aplus(real r2, real xr, real yr, real zr, real mscale, 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 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 f, real aewald, real &__restrict__ e, real &__restrict__ poti, real &__restrict__ potk, PairMPoleGrad &__restrict__ pgrad)
Definition: pairmpoleaplus.h:9
Generated by doxygen 1.9.5