Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
tortor.h
1#pragma once
2#include "math/const.h"
3#include "math/libfunc.h"
4#include "seq/add.h"
5#include "seq/seq.h"
6#include <tinker/detail/ktrtor.hh>
7
8namespace tinker {
9// see also bicubic.f
10#define BCUCOF__ \
11 real c[4][4]; \
12 real d1 = x1u - x1l; \
13 real d2 = x2u - x2l; \
14 real d12 = d1 * d2; \
15 real y1[4], y2[4], y12[4]; \
16 \
17 y1[0] = d1 * y1i[0]; \
18 y1[1] = d1 * y1i[1]; \
19 y1[2] = d1 * y1i[2]; \
20 y1[3] = d1 * y1i[3]; \
21 y2[0] = d2 * y2i[0]; \
22 y2[1] = d2 * y2i[1]; \
23 y2[2] = d2 * y2i[2]; \
24 y2[3] = d2 * y2i[3]; \
25 y12[0] = d12 * y12i[0]; \
26 y12[1] = d12 * y12i[1]; \
27 y12[2] = d12 * y12i[2]; \
28 y12[3] = d12 * y12i[3]; \
29 \
30 c[0][0] = y[0]; \
31 c[0][1] = y2[0]; \
32 c[0][2] = 3 * (y[3] - y[0]) - (2 * y2[0] + y2[3]); \
33 c[0][3] = 2 * (y[0] - y[3]) + y2[0] + y2[3]; \
34 c[1][0] = y1[0]; \
35 c[1][1] = y12[0]; \
36 c[1][2] = 3 * (y1[3] - y1[0]) - (2 * y12[0] + y12[3]); \
37 c[1][3] = 2 * (y1[0] - y1[3]) + y12[0] + y12[3]; \
38 c[2][0] = 3 * (y[1] - y[0]) - (2 * y1[0] + y1[1]); \
39 c[2][1] = 3 * (y2[1] - y2[0]) - (2 * y12[0] + y12[1]); \
40 c[2][2] = 9 * (y[0] - y[1] + y[2] - y[3]) + 6 * y1[0] + 3 * y1[1] \
41 - 3 * y1[2] - 6 * y1[3] + 6 * y2[0] - 6 * y2[1] - 3 * y2[2] + 3 * y2[3] \
42 + 4 * y12[0] + 2 * y12[1] + y12[2] + 2 * y12[3]; \
43 c[2][3] = 6 * (y[1] - y[0] + y[3] - y[2]) + -4 * y1[0] - 2 * y1[1] \
44 + 2 * y1[2] + 4 * y1[3] + -3 * y2[0] + 3 * y2[1] + 3 * y2[2] - 3 * y2[3] \
45 - 2 * y12[0] - y12[1] - y12[2] - 2 * y12[3]; \
46 c[3][0] = 2 * (y[0] - y[1]) + y1[0] + y1[1]; \
47 c[3][1] = 2 * (y2[0] - y2[1]) + y12[0] + y12[1]; \
48 c[3][2] = 6 * (y[1] - y[0] + y[3] - y[2]) \
49 + 3 * (y1[2] + y1[3] - y1[0] - y1[1]) \
50 + 2 * (2 * (y2[1] - y2[0]) + y2[2] - y2[3]) + -2 * (y12[0] + y12[1]) \
51 - y12[2] - y12[3]; \
52 c[3][3] = 4 * (y[0] - y[1] + y[2] - y[3]) \
53 + 2 * (y1[0] + y1[1] - y1[2] - y1[3]) \
54 + 2 * (y2[0] - y2[1] - y2[2] + y2[3]) + y12[0] + y12[1] + y12[2] \
55 + y12[3]; \
56 \
57 real t = (x1 - x1l) * REAL_RECIP(x1u - x1l); \
58 real u = (x2 - x2l) * REAL_RECIP(x2u - x2l)
59
60#define BCUCOF_ANSY__ \
61 real ay = ((c[3][3] * u + c[3][2]) * u + c[3][1]) * u + c[3][0]; \
62 ay = t * ay + ((c[2][3] * u + c[2][2]) * u + c[2][1]) * u + c[2][0]; \
63 ay = t * ay + ((c[1][3] * u + c[1][2]) * u + c[1][1]) * u + c[1][0]; \
64 ay = t * ay + ((c[0][3] * u + c[0][2]) * u + c[0][1]) * u + c[0][0]
65
66SEQ_ROUTINE
67static void bcuint0(const real (&restrict y)[4], const real (&restrict y1i)[4],
68 const real (&restrict y2i)[4], const real (&restrict y12i)[4], real x1l,
69 real x1u, real x2l, real x2u, real x1, real x2, real& restrict ansy)
70{
71 BCUCOF__;
72 BCUCOF_ANSY__;
73 ansy = ay;
74}
75
76SEQ_ROUTINE
77static void bcuint1(const real (&restrict y)[4], const real (&restrict y1i)[4],
78 const real (&restrict y2i)[4], const real (&restrict y12i)[4], real x1l,
79 real x1u, real x2l, real x2u, real x1, real x2, real& restrict ansy,
80 real& restrict ansy1, real& restrict ansy2)
81{
82 BCUCOF__;
83 BCUCOF_ANSY__;
84 ansy = ay;
85
86 real day1 = (3 * c[3][3] * t + 2 * c[2][3]) * t + c[1][3];
87 day1 = u * day1 + (3 * c[3][2] * t + 2 * c[2][2]) * t + c[1][2];
88 day1 = u * day1 + (3 * c[3][1] * t + 2 * c[2][1]) * t + c[1][1];
89 day1 = u * day1 + (3 * c[3][0] * t + 2 * c[2][0]) * t + c[1][0];
90 day1 *= REAL_RECIP(x1u - x1l);
91 ansy1 = day1;
92
93 real day2 = (3 * c[3][3] * u + 2 * c[3][2]) * u + c[3][1];
94 day2 = t * day2 + (3 * c[2][3] * u + 2 * c[2][2]) * u + c[2][1];
95 day2 = t * day2 + (3 * c[1][3] * u + 2 * c[1][2]) * u + c[1][1];
96 day2 = t * day2 + (3 * c[0][3] * u + 2 * c[0][2]) * u + c[0][1];
97 day2 *= REAL_RECIP(x2u - x2l);
98 ansy2 = day2;
99}
100#undef BCUCOF__
101#undef BCUCOF_ANSY__
102
103#pragma acc routine seq
104template <class Ver>
105SEQ_CUDA
106void dk_tortor(real& restrict e, real& restrict vxx, real& restrict vyx,
107 real& restrict vzx, real& restrict vyy, real& restrict vzy,
108 real& restrict vzz,
109
110 grad_prec* restrict dettx, grad_prec* restrict detty,
111 grad_prec* restrict dettz,
112
113 real ttorunit, int itortor, const int (*restrict itt)[3],
114 const int (*restrict ibitor)[5], const int* restrict chkttor_ia_,
115
116 const int* restrict tnx, const int* restrict tny,
117 const real (*restrict ttx)[ktrtor::maxtgrd],
118 const real (*restrict tty)[ktrtor::maxtgrd],
119 const real (*restrict tbf)[ktrtor::maxtgrd2],
120 const real (*restrict tbx)[ktrtor::maxtgrd2],
121 const real (*restrict tby)[ktrtor::maxtgrd2],
122 const real (*restrict tbxy)[ktrtor::maxtgrd2],
123
124 const real* restrict x, const real* restrict y, const real* restrict z)
125{
126 constexpr bool do_e = Ver::e;
127 constexpr bool do_g = Ver::g;
128 constexpr bool do_v = Ver::v;
129 if CONSTEXPR (do_e) e = 0;
130 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
131
132 real ftt[4], ft12[4], ft1[4], ft2[4];
133
134 int i = itt[itortor][0];
135 int k = itt[itortor][1]; // parameter indice
136 int ia, ib, ic, id, ie;
137 // if (itt(3,itortor) .eq. 1) then
138 // else
139 // end if
140 if (itt[itortor][2] == 0) {
141 ia = ibitor[i][0];
142 ib = ibitor[i][1];
143 ic = ibitor[i][2];
144 id = ibitor[i][3];
145 ie = ibitor[i][4];
146 } else {
147 ia = ibitor[i][4];
148 ib = ibitor[i][3];
149 ic = ibitor[i][2];
150 id = ibitor[i][1];
151 ie = ibitor[i][0];
152 }
153
154 // compute the values of the torsional angles
155
156 real xia = x[ia];
157 real yia = y[ia];
158 real zia = z[ia];
159 real xib = x[ib];
160 real yib = y[ib];
161 real zib = z[ib];
162 real xic = x[ic];
163 real yic = y[ic];
164 real zic = z[ic];
165 real xid = x[id];
166 real yid = y[id];
167 real zid = z[id];
168 real xie = x[ie];
169 real yie = y[ie];
170 real zie = z[ie];
171 real xba = xib - xia;
172 real yba = yib - yia;
173 real zba = zib - zia;
174 real xcb = xic - xib;
175 real ycb = yic - yib;
176 real zcb = zic - zib;
177 real xdc = xid - xic;
178 real ydc = yid - yic;
179 real zdc = zid - zic;
180 real xed = xie - xid;
181 real yed = yie - yid;
182 real zed = zie - zid;
183
184 real xt = yba * zcb - ycb * zba;
185 real yt = zba * xcb - zcb * xba;
186 real zt = xba * ycb - xcb * yba;
187 real xu = ycb * zdc - ydc * zcb;
188 real yu = zcb * xdc - zdc * xcb;
189 real zu = xcb * ydc - xdc * ycb;
190 real rt2 = xt * xt + yt * yt + zt * zt;
191 real ru2 = xu * xu + yu * yu + zu * zu;
192 real rtru = REAL_SQRT(rt2 * ru2);
193 real xv = ydc * zed - yed * zdc;
194 real yv = zdc * xed - zed * xdc;
195 real zv = xdc * yed - xed * ydc;
196 real rv2 = xv * xv + yv * yv + zv * zv;
197 real rurv = REAL_SQRT(ru2 * rv2);
198
199 if (rtru != 0 && rurv != 0) {
200 real sign;
201
202 real rcb = REAL_SQRT(xcb * xcb + ycb * ycb + zcb * zcb);
203 real cosine1 = (xt * xu + yt * yu + zt * zu) * REAL_RECIP(rtru);
204 cosine1 = REAL_MIN((real)1, REAL_MAX((real)-1, cosine1));
205 sign = xba * xu + yba * yu + zba * zu;
206 real angle1 = (sign < 0 ? -radian : radian) * REAL_ACOS(cosine1);
207 real value1 = angle1;
208
209 real rdc = REAL_SQRT(xdc * xdc + ydc * ydc + zdc * zdc);
210 real cosine2 = (xu * xv + yu * yv + zu * zv) * REAL_RECIP(rurv);
211 cosine2 = REAL_MIN((real)1, REAL_MAX((real)-1, cosine2));
212 sign = xcb * xv + ycb * yv + zcb * zv;
213 real angle2 = (sign < 0 ? -radian : radian) * REAL_ACOS(cosine2);
214 real value2 = angle2;
215
216 // check for inverted chirality at the central atom
217
218 const int chk_ia = chkttor_ia_[itortor];
219 real vol = 1;
220 if (chk_ia >= 0) {
221 // if chiral
222 real xac = x[chk_ia] - x[ic];
223 real yac = y[chk_ia] - y[ic];
224 real zac = z[chk_ia] - z[ic];
225 real c1 = -ycb * zdc + zcb * ydc;
226 real c2 = ydc * zac - zdc * yac;
227 real c3 = -yac * zcb + zac * ycb;
228 vol = xac * c1 - xcb * c2 + xdc * c3;
229 }
230 /*
231 * | non-chiral | chiral
232 * vol < 0 | sign = 1 | sign = -1; flip values
233 * vol >= 0 | sign = 1 | sign = 1
234 */
235 sign = (vol < 0 ? -1 : 1);
236 value1 = (vol < 0 ? -value1 : value1);
237 value2 = (vol < 0 ? -value2 : value2);
238 // angles must be within the range [-180, +180); +180 is not allowed
239 if (value1 < -180) value1 += 360;
240 if (value1 >= 180) value1 -= 360;
241 if (value2 < -180) value2 += 360;
242 if (value2 >= 180) value2 -= 360;
243
244 // use bicubic interpolation to compute spline values
245
246 // e.g., -180, -165, -150, ..., 165, 180, tnx[k] = 25
247 // xlo = floor((value+180) / (360/(25-1)))
248 int tnxk = tnx[k];
249 int xlo = REAL_FLOOR((value1 + 180) * (tnxk - 1) * REAL_RECIP(360));
250 int ylo = REAL_FLOOR((value2 + 180) * (tny[k] - 1) * REAL_RECIP(360));
251 real x1l = ttx[k][xlo];
252 real x1u = ttx[k][xlo + 1];
253 real y1l = tty[k][ylo];
254 real y1u = tty[k][ylo + 1];
255
256 int pos2 = (ylo + 1) * tnxk + xlo;
257 int pos1 = pos2 - tnxk;
258 ftt[0] = tbf[k][pos1];
259 ftt[1] = tbf[k][pos1 + 1];
260 ftt[2] = tbf[k][pos2 + 1];
261 ftt[3] = tbf[k][pos2];
262 ft1[0] = tbx[k][pos1];
263 ft1[1] = tbx[k][pos1 + 1];
264 ft1[2] = tbx[k][pos2 + 1];
265 ft1[3] = tbx[k][pos2];
266 ft2[0] = tby[k][pos1];
267 ft2[1] = tby[k][pos1 + 1];
268 ft2[2] = tby[k][pos2 + 1];
269 ft2[3] = tby[k][pos2];
270 ft12[0] = tbxy[k][pos1];
271 ft12[1] = tbxy[k][pos1 + 1];
272 ft12[2] = tbxy[k][pos2 + 1];
273 ft12[3] = tbxy[k][pos2];
274
275 real dedang1, dedang2;
276 if CONSTEXPR (do_g) {
277 bcuint1(ftt, ft1, ft2, ft12, x1l, x1u, y1l, y1u, value1, value2, e,
278 dedang1, dedang2);
279 } else {
280 bcuint0(ftt, ft1, ft2, ft12, x1l, x1u, y1l, y1u, value1, value2, e);
281 }
282
283 if CONSTEXPR (do_e) { e *= ttorunit; }
284
285 if CONSTEXPR (do_g) {
286 dedang1 *= (sign * ttorunit * radian);
287 dedang2 *= (sign * ttorunit * radian);
288
289 // chain rule terms for first angle derivative components
290
291 real xca = xic - xia;
292 real yca = yic - yia;
293 real zca = zic - zia;
294 real xdb = xid - xib;
295 real ydb = yid - yib;
296 real zdb = zid - zib;
297
298 real rt2cb_inv = REAL_RECIP(rt2 * rcb);
299 real ru2cb_inv = REAL_RECIP(ru2 * rcb);
300 real dedxt = dedang1 * (yt * zcb - ycb * zt) * rt2cb_inv;
301 real dedyt = dedang1 * (zt * xcb - zcb * xt) * rt2cb_inv;
302 real dedzt = dedang1 * (xt * ycb - xcb * yt) * rt2cb_inv;
303 real dedxu = -dedang1 * (yu * zcb - ycb * zu) * ru2cb_inv;
304 real dedyu = -dedang1 * (zu * xcb - zcb * xu) * ru2cb_inv;
305 real dedzu = -dedang1 * (xu * ycb - xcb * yu) * ru2cb_inv;
306
307 // compute first derivative components for first angle
308
309 real dedxia = zcb * dedyt - ycb * dedzt;
310 real dedyia = xcb * dedzt - zcb * dedxt;
311 real dedzia = ycb * dedxt - xcb * dedyt;
312 real dedxib = yca * dedzt - zca * dedyt + zdc * dedyu - ydc * dedzu;
313 real dedyib = zca * dedxt - xca * dedzt + xdc * dedzu - zdc * dedxu;
314 real dedzib = xca * dedyt - yca * dedxt + ydc * dedxu - xdc * dedyu;
315 real dedxic = zba * dedyt - yba * dedzt + ydb * dedzu - zdb * dedyu;
316 real dedyic = xba * dedzt - zba * dedxt + zdb * dedxu - xdb * dedzu;
317 real dedzic = yba * dedxt - xba * dedyt + xdb * dedyu - ydb * dedxu;
318 real dedxid = zcb * dedyu - ycb * dedzu;
319 real dedyid = xcb * dedzu - zcb * dedxu;
320 real dedzid = ycb * dedxu - xcb * dedyu;
321
322 // chain rule terms for second angle derivative components
323
324 real xec = xie - xic;
325 real yec = yie - yic;
326 real zec = zie - zic;
327
328 real ru2dc_inv = REAL_RECIP(ru2 * rdc);
329 real rv2dc_inv = REAL_RECIP(rv2 * rdc);
330 real dedxu2 = dedang2 * (yu * zdc - ydc * zu) * ru2dc_inv;
331 real dedyu2 = dedang2 * (zu * xdc - zdc * xu) * ru2dc_inv;
332 real dedzu2 = dedang2 * (xu * ydc - xdc * yu) * ru2dc_inv;
333 real dedxv2 = -dedang2 * (yv * zdc - ydc * zv) * rv2dc_inv;
334 real dedyv2 = -dedang2 * (zv * xdc - zdc * xv) * rv2dc_inv;
335 real dedzv2 = -dedang2 * (xv * ydc - xdc * yv) * rv2dc_inv;
336
337 // compute first derivative components for second angle
338
339 real dedxib2 = zdc * dedyu2 - ydc * dedzu2;
340 real dedyib2 = xdc * dedzu2 - zdc * dedxu2;
341 real dedzib2 = ydc * dedxu2 - xdc * dedyu2;
342 real dedxic2 = ydb * dedzu2 - zdb * dedyu2 + zed * dedyv2
343 - yed * dedzv2;
344 real dedyic2 = zdb * dedxu2 - xdb * dedzu2 + xed * dedzv2
345 - zed * dedxv2;
346 real dedzic2 = xdb * dedyu2 - ydb * dedxu2 + yed * dedxv2
347 - xed * dedyv2;
348 real dedxid2 = zcb * dedyu2 - ycb * dedzu2 + yec * dedzv2
349 - zec * dedyv2;
350 real dedyid2 = xcb * dedzu2 - zcb * dedxu2 + zec * dedxv2
351 - xec * dedzv2;
352 real dedzid2 = ycb * dedxu2 - xcb * dedyu2 + xec * dedyv2
353 - yec * dedxv2;
354 real dedxie2 = zdc * dedyv2 - ydc * dedzv2;
355 real dedyie2 = xdc * dedzv2 - zdc * dedxv2;
356 real dedzie2 = ydc * dedxv2 - xdc * dedyv2;
357
358 atomic_add(dedxia, dettx, ia);
359 atomic_add(dedyia, detty, ia);
360 atomic_add(dedzia, dettz, ia);
361 atomic_add(dedxib + dedxib2, dettx, ib);
362 atomic_add(dedyib + dedyib2, detty, ib);
363 atomic_add(dedzib + dedzib2, dettz, ib);
364 atomic_add(dedxic + dedxic2, dettx, ic);
365 atomic_add(dedyic + dedyic2, detty, ic);
366 atomic_add(dedzic + dedzic2, dettz, ic);
367 atomic_add(dedxid + dedxid2, dettx, id);
368 atomic_add(dedyid + dedyid2, detty, id);
369 atomic_add(dedzid + dedzid2, dettz, id);
370 atomic_add(dedxie2, dettx, ie);
371 atomic_add(dedyie2, detty, ie);
372 atomic_add(dedzie2, dettz, ie);
373
374 if CONSTEXPR (do_v) {
375 vxx = xcb * (dedxic + dedxid) - xba * dedxia + xdc * dedxid;
376 vyx = ycb * (dedxic + dedxid) - yba * dedxia + ydc * dedxid;
377 vzx = zcb * (dedxic + dedxid) - zba * dedxia + zdc * dedxid;
378 vyy = ycb * (dedyic + dedyid) - yba * dedyia + ydc * dedyid;
379 vzy = zcb * (dedyic + dedyid) - zba * dedyia + zdc * dedyid;
380 vzz = zcb * (dedzic + dedzid) - zba * dedzia + zdc * dedzid;
381 real vxx2 = xdc * (dedxid2 + dedxie2) - xcb * dedxib2
382 + xed * dedxie2;
383 real vyx2 = ydc * (dedxid2 + dedxie2) - ycb * dedxib2
384 + yed * dedxie2;
385 real vzx2 = zdc * (dedxid2 + dedxie2) - zcb * dedxib2
386 + zed * dedxie2;
387 real vyy2 = ydc * (dedyid2 + dedyie2) - ycb * dedyib2
388 + yed * dedyie2;
389 real vzy2 = zdc * (dedyid2 + dedyie2) - zcb * dedyib2
390 + zed * dedyie2;
391 real vzz2 = zdc * (dedzid2 + dedzie2) - zcb * dedzib2
392 + zed * dedzie2;
393
394 vxx += vxx2, vyx += vyx2, vzx += vzx2, vyy += vyy2, vzy += vzy2,
395 vzz += vzz2;
396 }
397 }
398 }
399}
400}
tinker::atomic_add
void atomic_add(T value, T *buffer, size_t offset=0)
Definition: acc/adddef.h:14
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::z
real * z
Current coordinates used in energy evaluation and neighbor lists.
tinker::x
real * x
Number of the trajectory frames.
tinker::y
real * y
Current coordinates used in energy evaluation and neighbor lists.
tinker::radian
constexpr real radian
Definition: const.h:14
tinker::real
float real
Definition: precision.h:80
tinker::grad_prec
fixed grad_prec
Definition: precision.h:103
tinker::ibitor
int(* ibitor)[5]
tinker::tby
real(* tby)[ktrtor::maxtgrd2]
tinker::tnx
int * tnx
tinker::dettx
grad_prec * dettx
tinker::itt
int(* itt)[3]
tinker::tbxy
real(* tbxy)[ktrtor::maxtgrd2]
tinker::ttorunit
real ttorunit
tinker::tny
int * tny
tinker::tty
real(* tty)[ktrtor::maxtgrd]
tinker::ttx
real(* ttx)[ktrtor::maxtgrd]
tinker::detty
grad_prec * detty
tinker::tbx
real(* tbx)[ktrtor::maxtgrd2]
tinker::chkttor_ia_
int * chkttor_ia_
tinker::dettz
grad_prec * dettz
tinker::tbf
real(* tbf)[ktrtor::maxtgrd2]
tinker
Definition: testrt.h:9
tinker::dk_tortor
__device__ void dk_tortor(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ dettx, grad_prec *__restrict__ detty, grad_prec *__restrict__ dettz, real ttorunit, int itortor, const int(*__restrict__ itt)[3], const int(*__restrict__ ibitor)[5], const int *__restrict__ chkttor_ia_, const int *__restrict__ tnx, const int *__restrict__ tny, const real(*__restrict__ ttx)[ktrtor::maxtgrd], const real(*__restrict__ tty)[ktrtor::maxtgrd], const real(*__restrict__ tbf)[ktrtor::maxtgrd2], const real(*__restrict__ tbx)[ktrtor::maxtgrd2], const real(*__restrict__ tby)[ktrtor::maxtgrd2], const real(*__restrict__ tbxy)[ktrtor::maxtgrd2], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: tortor.h:106
Generated by doxygen 1.9.5