Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
angtor.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
7namespace tinker {
8#pragma acc routine seq
9template <class Ver>
10SEQ_CUDA
11void dk_angtor(real& restrict e, real& restrict vxx, real& restrict vyx,
12 real& restrict vzx, real& restrict vyy, real& restrict vzy,
13 real& restrict vzz,
14
15 grad_prec* restrict deatx, grad_prec* restrict deaty,
16 grad_prec* restrict deatz,
17
18 real atorunit, int iangtor, const int (*restrict iat)[3],
19 const real (*restrict kant)[6],
20
21 const real* restrict anat, const int (*restrict itors)[4],
22 const real (*restrict tors1)[4], const real (*restrict tors2)[4],
23 const real (*restrict tors3)[4],
24
25 const real* restrict x, const real* restrict y, const real* restrict z)
26{
27 constexpr bool do_e = Ver::e;
28 constexpr bool do_g = Ver::g;
29 constexpr bool do_v = Ver::v;
30 if CONSTEXPR (do_e) e = 0;
31 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
32
33 const int i = iat[iangtor][0];
34 const int ia = itors[i][0];
35 const int ib = itors[i][1];
36 const int ic = itors[i][2];
37 const int id = itors[i][3];
38
39 real xia = x[ia];
40 real yia = y[ia];
41 real zia = z[ia];
42 real xib = x[ib];
43 real yib = y[ib];
44 real zib = z[ib];
45 real xic = x[ic];
46 real yic = y[ic];
47 real zic = z[ic];
48 real xid = x[id];
49 real yid = y[id];
50 real zid = z[id];
51 real xba = xib - xia;
52 real yba = yib - yia;
53 real zba = zib - zia;
54 real xcb = xic - xib;
55 real ycb = yic - yib;
56 real zcb = zic - zib;
57 real xdc = xid - xic;
58 real ydc = yid - yic;
59 real zdc = zid - zic;
60
61 real rba2 = xba * xba + yba * yba + zba * zba;
62 real rcb2 = xcb * xcb + ycb * ycb + zcb * zcb;
63 real rdc2 = xdc * xdc + ydc * ydc + zdc * zdc;
64 real rmin = REAL_MIN(rba2, rcb2);
65 rmin = REAL_MIN(rmin, rdc2);
66 if (rmin == 0) return;
67
68 real xt = yba * zcb - ycb * zba;
69 real yt = zba * xcb - zcb * xba;
70 real zt = xba * ycb - xcb * yba;
71 real xu = ycb * zdc - ydc * zcb;
72 real yu = zcb * xdc - zdc * xcb;
73 real zu = xcb * ydc - xdc * ycb;
74 real xtu = yt * zu - yu * zt;
75 real ytu = zt * xu - zu * xt;
76 real ztu = xt * yu - xu * yt;
77 real rt2 = xt * xt + yt * yt + zt * zt;
78 rt2 = REAL_MAX(rt2, (real)0.000001);
79 real ru2 = xu * xu + yu * yu + zu * zu;
80 ru2 = REAL_MAX(ru2, (real)0.000001);
81 real rtru = REAL_SQRT(rt2 * ru2);
82 real rcb = REAL_SQRT(rcb2);
83
84 real xca = xic - xia;
85 real yca = yic - yia;
86 real zca = zic - zia;
87 real xdb = xid - xib;
88 real ydb = yid - yib;
89 real zdb = zid - zib;
90
91 real cosine = (xt * xu + yt * yu + zt * zu) * REAL_RECIP(rtru);
92 real sine = (xcb * xtu + ycb * ytu + zcb * ztu) * REAL_RECIP(rcb * rtru);
93
94 real c1 = tors1[i][2];
95 real s1 = tors1[i][3];
96 real c2 = tors2[i][2];
97 real s2 = tors2[i][3];
98 real c3 = tors3[i][2];
99 real s3 = tors3[i][3];
100 real cosine2 = cosine * cosine - sine * sine;
101 real sine2 = 2 * cosine * sine;
102 real cosine3 = cosine * cosine2 - sine * sine2;
103 real sine3 = cosine * sine2 + sine * cosine2;
104 real phi1 = 1 + (cosine * c1 + sine * s1);
105 real phi2 = 1 + (cosine2 * c2 + sine2 * s2);
106 real phi3 = 1 + (cosine3 * c3 + sine3 * s3);
107
108 real dphi1, dphi2, dphi3;
109 if CONSTEXPR (do_g) {
110 dphi1 = cosine * s1 - sine * c1;
111 dphi2 = 2 * (cosine2 * s2 - sine2 * c2);
112 dphi3 = 3 * (cosine3 * s3 - sine3 * c3);
113 }
114
115 real e1, e2;
116 MAYBE_UNUSED real dedxia, dedyia, dedzia;
117 MAYBE_UNUSED real dedxib, dedyib, dedzib;
118 MAYBE_UNUSED real dedxic, dedyic, dedzic;
119 MAYBE_UNUSED real dedxid, dedyid, dedzid;
120
121 {
122 real v1 = kant[iangtor][0];
123 real v2 = kant[iangtor][1];
124 real v3 = kant[iangtor][2];
125 int k = iat[iangtor][1];
126 real dot = xba * xcb + yba * ycb + zba * zcb;
127 real cosang = -dot * REAL_RSQRT(rba2 * rcb2);
128 real angle = radian * REAL_ACOS(cosang);
129 real dt = angle - anat[k];
130 e1 = atorunit * dt * (v1 * phi1 + v2 * phi2 + v3 * phi3);
131 if CONSTEXPR (do_g) {
132 real dedphi = atorunit * dt * (v1 * dphi1 + v2 * dphi2 + v3 * dphi3);
133 real ddt = atorunit * radian * (v1 * phi1 + v2 * phi2 + v3 * phi3);
134 real dedxt = dedphi * (zcb * yt - ycb * zt) / (rt2 * rcb);
135 real dedyt = dedphi * (xcb * zt - zcb * xt) / (rt2 * rcb);
136 real dedzt = dedphi * (ycb * xt - xcb * yt) / (rt2 * rcb);
137 real dedxu = dedphi * (ycb * zu - zcb * yu) / (ru2 * rcb);
138 real dedyu = dedphi * (zcb * xu - xcb * zu) / (ru2 * rcb);
139 real dedzu = dedphi * (xcb * yu - ycb * xu) / (ru2 * rcb);
140
141 real terma = -ddt / (rba2 * REAL_SQRT(rt2));
142 real termc = ddt / (rcb2 * REAL_SQRT(rt2));
143 dedxia = terma * (zba * yt - yba * zt) + zcb * dedyt - ycb * dedzt;
144 dedyia = terma * (xba * zt - zba * xt) + xcb * dedzt - zcb * dedxt;
145 dedzia = terma * (yba * xt - xba * yt) + ycb * dedxt - xcb * dedyt;
146 dedxib = terma * (yba * zt - zba * yt) + termc * (zcb * yt - ycb * zt)
147 + yca * dedzt - zca * dedyt + zdc * dedyu - ydc * dedzu;
148 dedyib = terma * (zba * xt - xba * zt) + termc * (xcb * zt - zcb * xt)
149 + zca * dedxt - xca * dedzt + xdc * dedzu - zdc * dedxu;
150 dedzib = terma * (xba * yt - yba * xt) + termc * (ycb * xt - xcb * yt)
151 + xca * dedyt - yca * dedxt + ydc * dedxu - xdc * dedyu;
152 dedxic = termc * (ycb * zt - zcb * yt) + zba * dedyt - yba * dedzt
153 + ydb * dedzu - zdb * dedyu;
154 dedyic = termc * (zcb * xt - xcb * zt) + xba * dedzt - zba * dedxt
155 + zdb * dedxu - xdb * dedzu;
156 dedzic = termc * (xcb * yt - ycb * xt) + yba * dedxt - xba * dedyt
157 + xdb * dedyu - ydb * dedxu;
158 dedxid = zcb * dedyu - ycb * dedzu;
159 dedyid = xcb * dedzu - zcb * dedxu;
160 dedzid = ycb * dedxu - xcb * dedyu;
161 }
162 }
163
164 {
165 real v1 = kant[iangtor][3];
166 real v2 = kant[iangtor][4];
167 real v3 = kant[iangtor][5];
168 int k = iat[iangtor][2];
169 real dot = xcb * xdc + ycb * ydc + zcb * zdc;
170 real cosang = -dot * REAL_RSQRT(rcb2 * rdc2);
171 real angle = radian * REAL_ACOS(cosang);
172 real dt = angle - anat[k];
173 e2 = atorunit * dt * (v1 * phi1 + v2 * phi2 + v3 * phi3);
174 if CONSTEXPR (do_g) {
175 real dedphi = atorunit * dt * (v1 * dphi1 + v2 * dphi2 + v3 * dphi3);
176 real ddt = atorunit * radian * (v1 * phi1 + v2 * phi2 + v3 * phi3);
177 real dedxt = dedphi * (zcb * yt - ycb * zt) / (rt2 * rcb);
178 real dedyt = dedphi * (xcb * zt - zcb * xt) / (rt2 * rcb);
179 real dedzt = dedphi * (ycb * xt - xcb * yt) / (rt2 * rcb);
180 real dedxu = dedphi * (ycb * zu - zcb * yu) / (ru2 * rcb);
181 real dedyu = dedphi * (zcb * xu - xcb * zu) / (ru2 * rcb);
182 real dedzu = dedphi * (xcb * yu - ycb * xu) / (ru2 * rcb);
183
184 real termb = -ddt / (rcb2 * REAL_SQRT(ru2));
185 real termd = ddt / (rdc2 * REAL_SQRT(ru2));
186 dedxia += zcb * dedyt - ycb * dedzt;
187 dedyia += xcb * dedzt - zcb * dedxt;
188 dedzia += ycb * dedxt - xcb * dedyt;
189 dedxib += termb * (zcb * yu - ycb * zu) + yca * dedzt - zca * dedyt
190 + zdc * dedyu - ydc * dedzu;
191 dedyib += termb * (xcb * zu - zcb * xu) + zca * dedxt - xca * dedzt
192 + xdc * dedzu - zdc * dedxu;
193 dedzib += termb * (ycb * xu - xcb * yu) + xca * dedyt - yca * dedxt
194 + ydc * dedxu - xdc * dedyu;
195 dedxic += termb * (ycb * zu - zcb * yu) + termd * (zdc * yu - ydc * zu)
196 + zba * dedyt - yba * dedzt + ydb * dedzu - zdb * dedyu;
197 dedyic += termb * (zcb * xu - xcb * zu) + termd * (xdc * zu - zdc * xu)
198 + xba * dedzt - zba * dedxt + zdb * dedxu - xdb * dedzu;
199 dedzic += termb * (xcb * yu - ycb * xu) + termd * (ydc * xu - xdc * yu)
200 + yba * dedxt - xba * dedyt + xdb * dedyu - ydb * dedxu;
201 dedxid += termd * (ydc * zu - zdc * yu) + zcb * dedyu - ycb * dedzu;
202 dedyid += termd * (zdc * xu - xdc * zu) + xcb * dedzu - zcb * dedxu;
203 dedzid += termd * (xdc * yu - ydc * xu) + ycb * dedxu - xcb * dedyu;
204 }
205 }
206
207 if CONSTEXPR (do_e) { e = e1 + e2; }
208 if CONSTEXPR (do_g) {
209 atomic_add(dedxia, deatx, ia);
210 atomic_add(dedyia, deaty, ia);
211 atomic_add(dedzia, deatz, ia);
212 atomic_add(dedxib, deatx, ib);
213 atomic_add(dedyib, deaty, ib);
214 atomic_add(dedzib, deatz, ib);
215 atomic_add(dedxic, deatx, ic);
216 atomic_add(dedyic, deaty, ic);
217 atomic_add(dedzic, deatz, ic);
218 atomic_add(dedxid, deatx, id);
219 atomic_add(dedyid, deaty, id);
220 atomic_add(dedzid, deatz, id);
221 if CONSTEXPR (do_v) {
222 vxx = xcb * (dedxic + dedxid) - xba * dedxia + xdc * dedxid;
223 vyx = ycb * (dedxic + dedxid) - yba * dedxia + ydc * dedxid;
224 vzx = zcb * (dedxic + dedxid) - zba * dedxia + zdc * dedxid;
225 vyy = ycb * (dedyic + dedyid) - yba * dedyia + ydc * dedyid;
226 vzy = zcb * (dedyic + dedyid) - zba * dedyia + zdc * dedyid;
227 vzz = zcb * (dedzic + dedzid) - zba * dedzia + zdc * dedzid;
228 }
229 }
230}
231}
tinker::atomic_add
void atomic_add(T value, T *buffer, size_t offset=0)
Definition: acc/adddef.h:14
SEQ_CUDA
#define SEQ_CUDA
Definition: acc/seqdef.h:12
tinker::anat
real * anat
MAYBE_UNUSED
#define MAYBE_UNUSED
Definition: macro.h:75
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::kant
real(* kant)[6]
tinker::itors
int(* itors)[4]
tinker::deatz
grad_prec * deatz
tinker::tors2
real(* tors2)[4]
tinker::tors1
real(* tors1)[4]
tinker::atorunit
real atorunit
tinker::deaty
grad_prec * deaty
tinker::tors3
real(* tors3)[4]
tinker::iat
int(* iat)[3]
tinker::deatx
grad_prec * deatx
tinker
Definition: testrt.h:9
tinker::dk_angtor
__device__ void dk_angtor(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ deatx, grad_prec *__restrict__ deaty, grad_prec *__restrict__ deatz, real atorunit, int iangtor, const int(*__restrict__ iat)[3], const real(*__restrict__ kant)[6], const real *__restrict__ anat, const int(*__restrict__ itors)[4], const real(*__restrict__ tors1)[4], const real(*__restrict__ tors2)[4], const real(*__restrict__ tors3)[4], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: angtor.h:11
Generated by doxygen 1.9.5