Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pair_mpole_chgpen.h
1#pragma once
2#include "seq/damp_hippo.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( //
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 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_pole_v2<11>(dmpik, dmpi, dmpk, r, alphai, alphak);
74 } else {
75 damp_pole_v2<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 rr1i = bn[0] - (1 - mscale * dmpi[0]) * rr1;
102 rr3i = bn[1] - (1 - mscale * dmpi[1]) * rr3;
103 rr5i = bn[2] - (1 - mscale * dmpi[2]) * rr5;
104 rr7i = bn[3] - (1 - mscale * dmpi[3]) * rr7;
105 rr1k = bn[0] - (1 - mscale * dmpk[0]) * rr1;
106 rr3k = bn[1] - (1 - mscale * dmpk[1]) * rr3;
107 rr5k = bn[2] - (1 - mscale * dmpk[2]) * rr5;
108 rr7k = bn[3] - (1 - mscale * dmpk[3]) * rr7;
109 rr1ik = bn[0] - (1 - mscale * dmpik[0]) * rr1;
110 rr3ik = bn[1] - (1 - mscale * dmpik[1]) * rr3;
111 rr5ik = bn[2] - (1 - mscale * dmpik[2]) * rr5;
112 rr7ik = bn[3] - (1 - mscale * dmpik[3]) * rr7;
113 rr9ik = bn[4] - (1 - mscale * dmpik[4]) * rr9;
114 if CONSTEXPR (do_g) rr11ik = bn[5] - (1 - mscale * dmpik[5]) * rr11;
115 rr1 = bn[0] - (1 - mscale) * rr1;
116 rr3 = bn[1] - (1 - mscale) * rr3;
117
118 if CONSTEXPR (do_e)
119 e = term1 * rr1 + term4ik * rr7ik + term5ik * rr9ik + term1i * rr1i
120 + term1k * rr1k + term1ik * rr1ik + term2i * rr3i + term2k * rr3k
121 + term2ik * rr3ik + term3i * rr5i + term3k * rr5k + term3ik * rr5ik;
122 // end if (do_e)
123
124 if CONSTEXPR (do_g) {
125 // gradient
126 real qixk = qixx * qkx + qixy * qky + qixz * qkz;
127 real qiyk = qixy * qkx + qiyy * qky + qiyz * qkz;
128 real qizk = qixz * qkx + qiyz * qky + qizz * qkz;
129 real qkxi = qkxx * qix + qkxy * qiy + qkxz * qiz;
130 real qkyi = qkxy * qix + qkyy * qiy + qkyz * qiz;
131 real qkzi = qkxz * qix + qkyz * qiy + qkzz * qiz;
132
133 real diqkx = dix * qkxx + diy * qkxy + diz * qkxz;
134 real diqky = dix * qkxy + diy * qkyy + diz * qkyz;
135 real diqkz = dix * qkxz + diy * qkyz + diz * qkzz;
136 real dkqix = dkx * qixx + dky * qixy + dkz * qixz;
137 real dkqiy = dkx * qixy + dky * qiyy + dkz * qiyz;
138 real dkqiz = dkx * qixz + dky * qiyz + dkz * qizz;
139
140 real de = term1 * rr3 + term4ik * rr9ik + term5ik * rr11ik + term1i * rr3i
141 + term1k * rr3k + term1ik * rr3ik + term2i * rr5i + term2k * rr5k
142 + term2ik * rr5ik + term3i * rr7i + term3k * rr7k + term3ik * rr7ik;
143
144 term1 = -corek * rr3i - valk * rr3ik + dkr * rr5ik - qkr * rr7ik;
145 real term2 = corei * rr3k + vali * rr3ik + dir * rr5ik + qir * rr7ik;
146 real term3 = 2 * rr5ik;
147 real term4 = -2
148 * (corek * rr5i + valk * rr5ik - dkr * rr7ik + qkr * rr9ik);
149 real term5 = -2
150 * (corei * rr5k + vali * rr5ik + dir * rr7ik + qir * rr9ik);
151 real term6 = 4 * rr7ik;
152
153 if CONSTEXPR (CFLX) {
154 real t1i = corek * rr1i + valk * rr1ik;
155 real t1k = corei * rr1k + vali * rr1ik;
156 real t2i = -dkr * rr3ik;
157 real t2k = dir * rr3ik;
158 real t3i = qkr * rr5ik;
159 real t3k = qir * rr5ik;
160 poti = t1i + t2i + t3i;
161 potk = t1k + t2k + t3k;
162 }
163
164 pgrad.frcx = de * xr + term1 * dix + term2 * dkx + term3 * (diqkx - dkqix)
165 + term4 * qix + term5 * qkx + term6 * (qixk + qkxi);
166 pgrad.frcy = de * yr + term1 * diy + term2 * dky + term3 * (diqky - dkqiy)
167 + term4 * qiy + term5 * qky + term6 * (qiyk + qkyi);
168 pgrad.frcz = de * zr + term1 * diz + term2 * dkz + term3 * (diqkz - dkqiz)
169 + term4 * qiz + term5 * qkz + term6 * (qizk + qkzi);
170
171 // torque
172 real dirx = diy * zr - diz * yr;
173 real diry = diz * xr - dix * zr;
174 real dirz = dix * yr - diy * xr;
175 real dkrx = dky * zr - dkz * yr;
176 real dkry = dkz * xr - dkx * zr;
177 real dkrz = dkx * yr - dky * xr;
178 real dikx = diy * dkz - diz * dky;
179 real diky = diz * dkx - dix * dkz;
180 real dikz = dix * dky - diy * dkx;
181
182 real qirx = qiz * yr - qiy * zr;
183 real qiry = qix * zr - qiz * xr;
184 real qirz = qiy * xr - qix * yr;
185 real qkrx = qkz * yr - qky * zr;
186 real qkry = qkx * zr - qkz * xr;
187 real qkrz = qky * xr - qkx * yr;
188 real qikx = qky * qiz - qkz * qiy;
189 real qiky = qkz * qix - qkx * qiz;
190 real qikz = qkx * qiy - qky * qix;
191
192 real qikrx = qizk * yr - qiyk * zr;
193 real qikry = qixk * zr - qizk * xr;
194 real qikrz = qiyk * xr - qixk * yr;
195 real qkirx = qkzi * yr - qkyi * zr;
196 real qkiry = qkxi * zr - qkzi * xr;
197 real qkirz = qkyi * xr - qkxi * yr;
198
199 real diqkrx = diqkz * yr - diqky * zr;
200 real diqkry = diqkx * zr - diqkz * xr;
201 real diqkrz = diqky * xr - diqkx * yr;
202 real dkqirx = dkqiz * yr - dkqiy * zr;
203 real dkqiry = dkqix * zr - dkqiz * xr;
204 real dkqirz = dkqiy * xr - dkqix * yr;
205
206 real dqikx = diy * qkz - diz * qky + dky * qiz - dkz * qiy
207 - 2
208 * (qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy
209 - qiyz * qkyy - qizz * qkyz);
210 real dqiky = diz * qkx - dix * qkz + dkz * qix - dkx * qiz
211 - 2
212 * (qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz
213 - qixy * qkyz - qixz * qkzz);
214 real dqikz = dix * qky - diy * qkx + dkx * qiy - dky * qix
215 - 2
216 * (qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx
217 - qiyy * qkxy - qiyz * qkxz);
218
219 pgrad.ttmi[0] = -rr3ik * dikx + term1 * dirx + term3 * (dqikx + dkqirx)
220 - term4 * qirx - term6 * (qikrx + qikx);
221 pgrad.ttmi[1] = -rr3ik * diky + term1 * diry + term3 * (dqiky + dkqiry)
222 - term4 * qiry - term6 * (qikry + qiky);
223 pgrad.ttmi[2] = -rr3ik * dikz + term1 * dirz + term3 * (dqikz + dkqirz)
224 - term4 * qirz - term6 * (qikrz + qikz);
225 pgrad.ttmk[0] = rr3ik * dikx + term2 * dkrx - term3 * (dqikx + diqkrx)
226 - term5 * qkrx - term6 * (qkirx - qikx);
227 pgrad.ttmk[1] = rr3ik * diky + term2 * dkry - term3 * (dqiky + diqkry)
228 - term5 * qkry - term6 * (qkiry - qiky);
229 pgrad.ttmk[2] = rr3ik * dikz + term2 * dkrz - term3 * (dqikz + diqkrz)
230 - term5 * qkrz - term6 * (qkirz - qikz);
231 } // end if (do_g)
232}
233}
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
__device__ void pair_mpole_chgpen(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: pair_mpole_chgpen.h:9
Generated by doxygen 1.9.5