Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
pair_repel.h
1#pragma once
2#include "math/realn.h"
3#include "math/switch.h"
4#include "seq/damp_hippo.h"
5#include "seq/pair_vlambda.h"
6
7namespace tinker {
8struct PairRepelGrad
9{
10 real frcx, frcy, frcz;
11 real ttqi[3];
12 real ttqk[3];
13};
14
15SEQ_ROUTINE
16inline void zero(PairRepelGrad& pgrad)
17{
18 pgrad.frcx = 0;
19 pgrad.frcy = 0;
20 pgrad.frcz = 0;
21 pgrad.ttqi[0] = 0;
22 pgrad.ttqi[1] = 0;
23 pgrad.ttqi[2] = 0;
24 pgrad.ttqk[0] = 0;
25 pgrad.ttqk[1] = 0;
26 pgrad.ttqk[2] = 0;
27}
28
29#pragma acc routine seq
30template <bool do_g, int SOFTCORE>
31SEQ_CUDA
32void pair_repel(real r2,
33 real rscale,
34 real vlambda,
35 real cut,
36 real off,
37 real xr,
38 real yr,
39 real zr,
40 real sizi,
41 real dmpi,
42 real vali,
43 real ci,
44 real dix,
45 real diy,
46 real diz,
47 real qixx,
48 real qixy,
49 real qixz,
50 real qiyy,
51 real qiyz,
52 real qizz,
53 real sizk,
54 real dmpk,
55 real valk,
56 real ck,
57 real dkx,
58 real dky,
59 real dkz,
60 real qkxx,
61 real qkxy,
62 real qkxz,
63 real qkyy,
64 real qkyz,
65 real qkzz,
66 real& restrict e,
67 PairRepelGrad& restrict pgrad)
68{
69 real cut2 = cut * cut;
70 real r = REAL_SQRT(r2);
71 real rInv = REAL_RECIP(r);
72 real rr2 = rInv * rInv;
73 real rr1 = rInv;
74 real rr3 = rr1 * rr2;
75 real rr5 = 3 * rr3 * rr2;
76 real rr7 = 5 * rr5 * rr2;
77 real rr9 = 7 * rr7 * rr2;
78 real rr11 = 0;
79 if CONSTEXPR (do_g)
80 rr11 = 9 * rr9 * rr2;
81
82 real3 dr = make_real3(xr, yr, zr);
83 real3 di = make_real3(dix, diy, diz);
84 real3 dk = make_real3(dkx, dky, dkz);
85
86 // get damping coefficients for the Pauli repulsion energy
87 real dmpik[6];
88
89 if CONSTEXPR (do_g) {
90 damp_rep<11>(dmpik, r, rr1, r2, rr3, rr5, rr7, rr9, rr11, dmpi, dmpk);
91 } else {
92 damp_rep<9>(dmpik, r, rr1, r2, rr3, rr5, rr7, rr9, rr11, dmpi, dmpk);
93 }
94
95 // Fortran tinker indexing
96 // dmpik(1) == dmpik[0]
97 // dmpik(3) == dmpik[1]
98 // dmpik(5) == dmpik[2]
99 // dmpik(7) == dmpik[3]
100 // dmpik(9) == dmpik[4]
101 // dmpik(11) == dmpik[5]
102
103 // same as in emplar_cu
104 real dir = dot3(di, dr);
105 real3 qi_dr = matvec(qixx, qixy, qixz, qiyy, qiyz, qizz, dr);
106 real qir = dot3(dr, qi_dr);
107 real dkr = dot3(dk, dr);
108 real3 qk_dr = matvec(qkxx, qkxy, qkxz, qkyy, qkyz, qkzz, dr);
109 real qkr = dot3(dr, qk_dr);
110
111 real dik = dot3(di, dk);
112 real qik = dot3(qi_dr, qk_dr);
113 real diqk = dot3(di, qk_dr);
114 real dkqi = dot3(dk, qi_dr);
115 real qiqk = dot3(qixx, qiyy, qizz, qkxx, qkyy, qkzz) + 2 * dot3(qixy, qixz, qiyz, qkxy, qkxz, qkyz);
116
117 real term1 = vali * valk;
118 real term2 = valk * dir - vali * dkr + dik;
119 real term3 = vali * qkr + valk * qir - dir * dkr + 2 * (dkqi - diqk + qiqk);
120 real term4 = dir * qkr - dkr * qir - 4 * qik;
121 real term5 = qir * qkr;
122
123 real sizik = sizi * sizk * rscale;
124 real eterm = term1 * dmpik[0] + term2 * dmpik[1] + term3 * dmpik[2] + term4 * dmpik[3] + term5 * dmpik[4];
125
126 e = sizik * eterm * rInv;
127
128 // energy via soft core lambda scaling
129 real termsc, soft;
130 if CONSTEXPR (SOFTCORE) {
131 real vlambda3 = vlambda * vlambda * vlambda;
132 real vlambda4 = vlambda3 * vlambda;
133 real vlambda5 = vlambda4 * vlambda;
134 termsc = vlambda3 - vlambda4 + r2;
135 soft = vlambda5 * r / REAL_SQRT(termsc);
136 e *= soft;
137 }
138 // gradient
139 if CONSTEXPR (do_g) {
140 real de = term1 * dmpik[1] + term2 * dmpik[2] + term3 * dmpik[3] + term4 * dmpik[4] + term5 * dmpik[5];
141 term1 = -valk * dmpik[1] + dkr * dmpik[2] - qkr * dmpik[3];
142 term2 = vali * dmpik[1] + dir * dmpik[2] + qir * dmpik[3];
143 term3 = 2 * dmpik[2];
144 term4 = 2 * (-valk * dmpik[2] + dkr * dmpik[3] - qkr * dmpik[4]);
145 term5 = 2 * (-vali * dmpik[2] - dir * dmpik[3] - qir * dmpik[4]);
146 real term6 = 4 * dmpik[3];
147
148 // few intermediates, same as in emplar
149 real qix = qi_dr.x;
150 real qiy = qi_dr.y;
151 real qiz = qi_dr.z;
152 real qkx = qk_dr.x;
153 real qky = qk_dr.y;
154 real qkz = qk_dr.z;
155
156 real qixk = qixx * qkx + qixy * qky + qixz * qkz; // |Qi Qk r>
157 real qiyk = qixy * qkx + qiyy * qky + qiyz * qkz;
158 real qizk = qixz * qkx + qiyz * qky + qizz * qkz;
159 real qkxi = qkxx * qix + qkxy * qiy + qkxz * qiz; // |Qk Qi r>
160 real qkyi = qkxy * qix + qkyy * qiy + qkyz * qiz;
161 real qkzi = qkxz * qix + qkyz * qiy + qkzz * qiz;
162
163 real diqkx = di.x * qkxx + di.y * qkxy + di.z * qkxz; // |Qk Di>
164 real diqky = di.x * qkxy + di.y * qkyy + di.z * qkyz;
165 real diqkz = di.x * qkxz + di.y * qkyz + di.z * qkzz;
166 real dkqix = dk.x * qixx + dk.y * qixy + dk.z * qixz; // |Qi Dk>
167 real dkqiy = dk.x * qixy + dk.y * qiyy + dk.z * qiyz;
168 real dkqiz = dk.x * qixz + dk.y * qiyz + dk.z * qizz;
169
170 real3 frc0;
171 frc0.x = de * dr.x + term1 * di.x + term2 * dk.x + term3 * (diqkx - dkqix) + term4 * qix + term5 * qkx + term6 * (qixk + qkxi);
172 frc0.y = de * dr.y + term1 * di.y + term2 * dk.y + term3 * (diqky - dkqiy) + term4 * qiy + term5 * qky + term6 * (qiyk + qkyi);
173 frc0.z = de * dr.z + term1 * di.z + term2 * dk.z + term3 * (diqkz - dkqiz) + term4 * qiz + term5 * qkz + term6 * (qizk + qkzi);
174
175 pgrad.frcx = frc0.x * rr1 + eterm * rr3 * dr.x;
176 pgrad.frcy = frc0.y * rr1 + eterm * rr3 * dr.y;
177 pgrad.frcz = frc0.z * rr1 + eterm * rr3 * dr.z;
178
179 pgrad.frcx *= sizik;
180 pgrad.frcy *= sizik;
181 pgrad.frcz *= sizik;
182
183 // torque
184 real dirx = di.y * dr.z - di.z * dr.y; // Di x r
185 real diry = di.z * dr.x - di.x * dr.z;
186 real dirz = di.x * dr.y - di.y * dr.x;
187 real dkrx = dk.y * dr.z - dk.z * dr.y; // Dk x r
188 real dkry = dk.z * dr.x - dk.x * dr.z;
189 real dkrz = dk.x * dr.y - dk.y * dr.x;
190 real dikx = di.y * dk.z - di.z * dk.y; // Di x Dk
191 real diky = di.z * dk.x - di.x * dk.z;
192 real dikz = di.x * dk.y - di.y * dk.x;
193
194 real qirx = qiz * dr.y - qiy * dr.z; // r x (Qi r)
195 real qiry = qix * dr.z - qiz * dr.x;
196 real qirz = qiy * dr.x - qix * dr.y;
197 real qkrx = qkz * dr.y - qky * dr.z; // r x (Qk r)
198 real qkry = qkx * dr.z - qkz * dr.x;
199 real qkrz = qky * dr.x - qkx * dr.y;
200 real qikx = qky * qiz - qkz * qiy; // (Qk r) x (Qi r)
201 real qiky = qkz * qix - qkx * qiz;
202 real qikz = qkx * qiy - qky * qix;
203
204 real qikrx = qizk * dr.y - qiyk * dr.z;
205 real qikry = qixk * dr.z - qizk * dr.x;
206 real qikrz = qiyk * dr.x - qixk * dr.y;
207 real qkirx = qkzi * dr.y - qkyi * dr.z;
208 real qkiry = qkxi * dr.z - qkzi * dr.x;
209 real qkirz = qkyi * dr.x - qkxi * dr.y;
210
211 real diqkrx = diqkz * dr.y - diqky * dr.z;
212 real diqkry = diqkx * dr.z - diqkz * dr.x;
213 real diqkrz = diqky * dr.x - diqkx * dr.y;
214 real dkqirx = dkqiz * dr.y - dkqiy * dr.z;
215 real dkqiry = dkqix * dr.z - dkqiz * dr.x;
216 real dkqirz = dkqiy * dr.x - dkqix * dr.y;
217
218 real dqikx = di.y * qkz - di.z * qky + dk.y * qiz - dk.z * qiy - 2 * (qixy * qkxz + qiyy * qkyz + qiyz * qkzz - qixz * qkxy - qiyz * qkyy - qizz * qkyz);
219 real dqiky = di.z * qkx - di.x * qkz + dk.z * qix - dk.x * qiz - 2 * (qixz * qkxx + qiyz * qkxy + qizz * qkxz - qixx * qkxz - qixy * qkyz - qixz * qkzz);
220 real dqikz = di.x * qky - di.y * qkx + dk.x * qiy - dk.y * qix - 2 * (qixx * qkxy + qixy * qkyy + qixz * qkyz - qixy * qkxx - qiyy * qkxy - qiyz * qkxz);
221
222 real3 trq0;
223 trq0.x = -dmpik[1] * dikx + term1 * dirx + term3 * (dqikx + dkqirx) - term4 * qirx - term6 * (qikrx + qikx);
224 trq0.y = -dmpik[1] * diky + term1 * diry + term3 * (dqiky + dkqiry) - term4 * qiry - term6 * (qikry + qiky);
225 trq0.z = -dmpik[1] * dikz + term1 * dirz + term3 * (dqikz + dkqirz) - term4 * qirz - term6 * (qikrz + qikz);
226
227 pgrad.ttqi[0] += sizik * rr1 * trq0.x;
228 pgrad.ttqi[1] += sizik * rr1 * trq0.y;
229 pgrad.ttqi[2] += sizik * rr1 * trq0.z;
230
231 trq0.x = dmpik[1] * dikx + term2 * dkrx - term3 * (dqikx + diqkrx) - term5 * qkrx - term6 * (qkirx - qikx);
232 trq0.y = dmpik[1] * diky + term2 * dkry - term3 * (dqiky + diqkry) - term5 * qkry - term6 * (qkiry - qiky);
233 trq0.z = dmpik[1] * dikz + term2 * dkrz - term3 * (dqikz + diqkrz) - term5 * qkrz - term6 * (qkirz - qikz);
234
235 pgrad.ttqk[0] += sizik * rr1 * trq0.x;
236 pgrad.ttqk[1] += sizik * rr1 * trq0.y;
237 pgrad.ttqk[2] += sizik * rr1 * trq0.z;
238
239 // force via soft core lambda scaling
240 if CONSTEXPR (SOFTCORE) {
241 real dsoft = e * soft * (rr2 - 1 / termsc);
242 pgrad.frcx = pgrad.frcx * soft - dsoft * xr;
243 pgrad.frcy = pgrad.frcy * soft - dsoft * yr;
244 pgrad.frcz = pgrad.frcz * soft - dsoft * zr;
245
246 pgrad.ttqi[0] *= soft;
247 pgrad.ttqi[1] *= soft;
248 pgrad.ttqi[2] *= soft;
249 pgrad.ttqk[0] *= soft;
250 pgrad.ttqk[1] *= soft;
251 pgrad.ttqk[2] *= soft;
252 }
253 }
254
255 if (r2 > cut2) {
256 real taper, dtaper;
257 switchTaper5<do_g>(r, cut, off, taper, dtaper);
258 if CONSTEXPR (do_g) {
259 dtaper *= e * rr1;
260 pgrad.frcx = pgrad.frcx * taper - dtaper * xr;
261 pgrad.frcy = pgrad.frcy * taper - dtaper * yr;
262 pgrad.frcz = pgrad.frcz * taper - dtaper * zr;
263
264 pgrad.ttqi[0] *= taper;
265 pgrad.ttqi[1] *= taper;
266 pgrad.ttqi[2] *= taper;
267 pgrad.ttqk[0] *= taper;
268 pgrad.ttqk[1] *= taper;
269 pgrad.ttqk[2] *= taper;
270 }
271 e *= taper;
272 }
273}
274}
tinker::float3::x
float x
Definition: acc/realndef.h:49
tinker::float3::y
float y
Definition: acc/realndef.h:49
tinker::float3::z
float z
Definition: acc/realndef.h:49
SEQ_ROUTINE
#define SEQ_ROUTINE
Definition: acc/seqdef.h:7
SEQ_CUDA
#define SEQ_CUDA
Definition: acc/seqdef.h:12
tinker::float3
Definition: acc/realndef.h:48
restrict
#define restrict
Definition: macro.h:51
CONSTEXPR
#define CONSTEXPR
Definition: macro.h:61
tinker::dot3
__device__ real dot3(real3 a, real3 b)
Definition: realn.h:84
tinker::matvec
__device__ real3 matvec(real xx, real xy, real xz, real yy, real yz, real zz, real3 v)
Definition: realn.h:109
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_repel
__device__ void pair_repel(real r2, real rscale, real vlambda, real cut, real off, real xr, real yr, real zr, real sizi, real dmpi, real vali, real ci, real dix, real diy, real diz, real qixx, real qixy, real qixz, real qiyy, real qiyz, real qizz, real sizk, real dmpk, real valk, real ck, real dkx, real dky, real dkz, real qkxx, real qkxy, real qkxz, real qkyy, real qkyz, real qkzz, real &__restrict__ e, PairRepelGrad &__restrict__ pgrad)
Definition: pair_repel.h:32
tinker::PairRepelGrad
Definition: pair_repel.h:9
tinker::PairRepelGrad::ttqk
real ttqk[3]
Definition: pair_repel.h:12
tinker::PairRepelGrad::ttqi
real ttqi[3]
Definition: pair_repel.h:11
tinker::PairRepelGrad::frcy
real frcy
Definition: pair_repel.h:10
tinker::PairRepelGrad::frcz
real frcz
Definition: pair_repel.h:10
tinker::PairRepelGrad::frcx
real frcx
Definition: pair_repel.h:10
Generated by doxygen 1.9.5