Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
angle.h
1#pragma once
2#include "ff/evalence.h"
3#include "math/const.h"
4#include "math/libfunc.h"
5#include "seq/add.h"
6#include "seq/seq.h"
7
8namespace tinker {
26#pragma acc routine seq
27template <class Ver>
28SEQ_CUDA
29void dk_angle(real& restrict e, real& restrict vxx, real& restrict vyx,
30 real& restrict vzx, real& restrict vyy, real& restrict vzy,
31 real& restrict vzz,
32
33 grad_prec* restrict deax, grad_prec* restrict deay, grad_prec* restrict deaz,
34
35 const Angle* restrict angtyp, real angunit, int i,
36 const int (*restrict iang)[4], const real* restrict anat,
37 const real* restrict ak, const real* restrict afld,
38
39 real cang, real qang, real pang, real sang,
40
41 const real* restrict x, const real* restrict y, const real* restrict z)
42{
43 constexpr bool do_e = Ver::e;
44 constexpr bool do_g = Ver::g;
45 constexpr bool do_v = Ver::v;
46 if CONSTEXPR (do_e) e = 0;
47 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
48
49 int ia = iang[i][0];
50 int ib = iang[i][1];
51 int ic = iang[i][2];
52 int id = iang[i][3];
53 real ideal = anat[i];
54 real force = ak[i];
55 Angle angtypi = angtyp[i];
56
57 real xia = x[ia];
58 real yia = y[ia];
59 real zia = z[ia];
60 real xib = x[ib];
61 real yib = y[ib];
62 real zib = z[ib];
63 real xic = x[ic];
64 real yic = y[ic];
65 real zic = z[ic];
66
67 if (angtypi != Angle::IN_PLANE) {
68 real xab = xia - xib;
69 real yab = yia - yib;
70 real zab = zia - zib;
71 real xcb = xic - xib;
72 real ycb = yic - yib;
73 real zcb = zic - zib;
74
75 real rab2 = xab * xab + yab * yab + zab * zab;
76 real rcb2 = xcb * xcb + ycb * ycb + zcb * zcb;
77
78 if (rab2 != 0 && rcb2 != 0) {
79 real xp = ycb * zab - zcb * yab;
80 real yp = zcb * xab - xcb * zab;
81 real zp = xcb * yab - ycb * xab;
82 real rp = REAL_SQRT(xp * xp + yp * yp + zp * zp);
83 rp = REAL_MAX(rp, (real)0.0001);
84 real dot = xab * xcb + yab * ycb + zab * zcb;
85 real cosine = dot * REAL_RSQRT(rab2 * rcb2);
86 cosine = REAL_MIN((real)1, REAL_MAX((real)-1, cosine));
87 real angle = radian * REAL_ACOS(cosine);
88
89 real deddt;
90 if (angtypi == Angle::HARMONIC) {
91 real dt = angle - ideal;
92 real dt2 = dt * dt;
93 real dt3 = dt2 * dt;
94 real dt4 = dt2 * dt2;
95 if CONSTEXPR (do_e)
96 e = angunit * force * dt2
97 * (1 + cang * dt + qang * dt2 + pang * dt3 + sang * dt4);
98 if CONSTEXPR (do_g)
99 deddt = angunit * force * dt * radian
100 * (2 + 3 * cang * dt + 4 * qang * dt2 + 5 * pang * dt3
101 + 6 * sang * dt4);
102 } else if (angtypi == Angle::LINEAR) {
103 real factor = 2 * angunit * radian * radian;
104 real sine = REAL_SQRT(1 - cosine * cosine);
105 if CONSTEXPR (do_e) e = factor * force * (1 + cosine);
106 if CONSTEXPR (do_g) deddt = -factor * force * sine;
107 } else if (angtypi == Angle::FOURIER) {
108 real fold = afld[i];
109 real factor = 2 * angunit * (radian / fold) * (radian / fold);
110 real dt = (fold * angle - ideal) * _1radian;
111 if CONSTEXPR (do_e) {
112 real cosine = REAL_COS(dt);
113 e = factor * force * (1 + cosine);
114 }
115 if CONSTEXPR (do_g) {
116 real sine = REAL_SIN(dt);
117 deddt = -factor * force * fold * sine;
118 }
119 }
120
121 if CONSTEXPR (do_g) {
122 real terma = -deddt * REAL_RECIP(rab2 * rp);
123 real termc = deddt * REAL_RECIP(rcb2 * rp);
124 real dedxia = terma * (yab * zp - zab * yp);
125 real dedyia = terma * (zab * xp - xab * zp);
126 real dedzia = terma * (xab * yp - yab * xp);
127 real dedxic = termc * (ycb * zp - zcb * yp);
128 real dedyic = termc * (zcb * xp - xcb * zp);
129 real dedzic = termc * (xcb * yp - ycb * xp);
130 real dedxib = -dedxia - dedxic;
131 real dedyib = -dedyia - dedyic;
132 real dedzib = -dedzia - dedzic;
133
134 atomic_add(dedxia, deax, ia);
135 atomic_add(dedyia, deay, ia);
136 atomic_add(dedzia, deaz, ia);
137 atomic_add(dedxib, deax, ib);
138 atomic_add(dedyib, deay, ib);
139 atomic_add(dedzib, deaz, ib);
140 atomic_add(dedxic, deax, ic);
141 atomic_add(dedyic, deay, ic);
142 atomic_add(dedzic, deaz, ic);
143
144 if CONSTEXPR (do_v) {
145 vxx = xab * dedxia + xcb * dedxic;
146 vyx = yab * dedxia + ycb * dedxic;
147 vzx = zab * dedxia + zcb * dedxic;
148 vyy = yab * dedyia + ycb * dedyic;
149 vzy = zab * dedyia + zcb * dedyic;
150 vzz = zab * dedzia + zcb * dedzic;
151 }
152 }
153 }
154 } else {
155 real xid = x[id];
156 real yid = y[id];
157 real zid = z[id];
158 real xad = xia - xid;
159 real yad = yia - yid;
160 real zad = zia - zid;
161 real xbd = xib - xid;
162 real ybd = yib - yid;
163 real zbd = zib - zid;
164 real xcd = xic - xid;
165 real ycd = yic - yid;
166 real zcd = zic - zid;
167 real xt = yad * zcd - zad * ycd;
168 real yt = zad * xcd - xad * zcd;
169 real zt = xad * ycd - yad * xcd;
170 real rt2 = xt * xt + yt * yt + zt * zt;
171 real delta = -(xt * xbd + yt * ybd + zt * zbd) * REAL_RECIP(rt2);
172 real xap = xia - xib - xt * delta;
173 real yap = yia - yib - yt * delta;
174 real zap = zia - zib - zt * delta;
175 real xcp = xic - xib - xt * delta;
176 real ycp = yic - yib - yt * delta;
177 real zcp = zic - zib - zt * delta;
178 real rap2 = xap * xap + yap * yap + zap * zap;
179 real rcp2 = xcp * xcp + ycp * ycp + zcp * zcp;
180 if (rap2 != 0 && rcp2 != 0) {
181 real dot = xap * xcp + yap * ycp + zap * zcp;
182 real cosine = dot * REAL_RSQRT(rap2 * rcp2);
183 cosine = REAL_MIN((real)1, REAL_MAX((real)-1, cosine));
184 real angle = radian * REAL_ACOS(cosine);
185 real dt = angle - ideal;
186 real dt2 = dt * dt;
187 real dt3 = dt2 * dt;
188 real dt4 = dt2 * dt2;
189
190 if CONSTEXPR (do_e) {
191 e = angunit * force * dt2
192 * (1 + cang * dt + qang * dt2 + pang * dt3 + sang * dt4);
193 }
194
195 if CONSTEXPR (do_g) {
196 real deddt = angunit * force * dt * radian
197 * (2 + 3 * cang * dt + 4 * qang * dt2 + 5 * pang * dt3
198 + 6 * sang * dt4);
199 real xm = ycp * zap - zcp * yap;
200 real ym = zcp * xap - xcp * zap;
201 real zm = xcp * yap - ycp * xap;
202 real rm = REAL_SQRT(xm * xm + ym * ym + zm * zm);
203 rm = REAL_MAX(rm, (real)0.0001);
204
205 // chain rule terms for first derivative components
206
207 real terma = -deddt * REAL_RECIP(rap2 * rm);
208 real termc = deddt * REAL_RECIP(rcp2 * rm);
209 real dedxia = terma * (yap * zm - zap * ym);
210 real dedyia = terma * (zap * xm - xap * zm);
211 real dedzia = terma * (xap * ym - yap * xm);
212 real dedxic = termc * (ycp * zm - zcp * ym);
213 real dedyic = termc * (zcp * xm - xcp * zm);
214 real dedzic = termc * (xcp * ym - ycp * xm);
215 real dedxip = -dedxia - dedxic;
216 real dedyip = -dedyia - dedyic;
217 real dedzip = -dedzia - dedzic;
218
219 // chain rule components for the projection of the central atom
220
221 real delta2, term;
222 delta2 = 2 * delta;
223 real ptrt2 = (dedxip * xt + dedyip * yt + dedzip * zt)
224 * REAL_RECIP(rt2);
225 term = (zcd * ybd - ycd * zbd) + delta2 * (yt * zcd - zt * ycd);
226 real dpdxia = delta * (ycd * dedzip - zcd * dedyip) + term * ptrt2;
227 term = (xcd * zbd - zcd * xbd) + delta2 * (zt * xcd - xt * zcd);
228 real dpdyia = delta * (zcd * dedxip - xcd * dedzip) + term * ptrt2;
229 term = (ycd * xbd - xcd * ybd) + delta2 * (xt * ycd - yt * xcd);
230 real dpdzia = delta * (xcd * dedyip - ycd * dedxip) + term * ptrt2;
231 term = (yad * zbd - zad * ybd) + delta2 * (zt * yad - yt * zad);
232 real dpdxic = delta * (zad * dedyip - yad * dedzip) + term * ptrt2;
233 term = (zad * xbd - xad * zbd) + delta2 * (xt * zad - zt * xad);
234 real dpdyic = delta * (xad * dedzip - zad * dedxip) + term * ptrt2;
235 term = (xad * ybd - yad * xbd) + delta2 * (yt * xad - xt * yad);
236 real dpdzic = delta * (yad * dedxip - xad * dedyip) + term * ptrt2;
237
238 // compute derivative components for this interaction
239
240 dedxia = dedxia + dpdxia;
241 dedyia = dedyia + dpdyia;
242 dedzia = dedzia + dpdzia;
243 real dedxib = dedxip;
244 real dedyib = dedyip;
245 real dedzib = dedzip;
246 dedxic = dedxic + dpdxic;
247 dedyic = dedyic + dpdyic;
248 dedzic = dedzic + dpdzic;
249 real dedxid = -dedxia - dedxib - dedxic;
250 real dedyid = -dedyia - dedyib - dedyic;
251 real dedzid = -dedzia - dedzib - dedzic;
252
253 atomic_add(dedxia, deax, ia);
254 atomic_add(dedyia, deay, ia);
255 atomic_add(dedzia, deaz, ia);
256 atomic_add(dedxib, deax, ib);
257 atomic_add(dedyib, deay, ib);
258 atomic_add(dedzib, deaz, ib);
259 atomic_add(dedxic, deax, ic);
260 atomic_add(dedyic, deay, ic);
261 atomic_add(dedzic, deaz, ic);
262 atomic_add(dedxid, deax, id);
263 atomic_add(dedyid, deay, id);
264 atomic_add(dedzid, deaz, id);
265
266 if CONSTEXPR (do_v) {
267 vxx = xad * dedxia + xbd * dedxib + xcd * dedxic;
268 vyx = yad * dedxia + ybd * dedxib + ycd * dedxic;
269 vzx = zad * dedxia + zbd * dedxib + zcd * dedxic;
270 vyy = yad * dedyia + ybd * dedyib + ycd * dedyic;
271 vzy = zad * dedyia + zbd * dedyib + zcd * dedyic;
272 vzz = zad * dedzia + zbd * dedzib + zcd * dedzic;
273 }
274 }
275 }
276 }
277}
278}
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
tinker::afld
real * afld
tinker::deay
grad_prec * deay
tinker::sang
real sang
tinker::angunit
real angunit
tinker::qang
real qang
tinker::angtyp
Angle * angtyp
tinker::iang
int(* iang)[4]
tinker::Angle
Angle
Definition: evalence.h:68
tinker::ak
real * ak
tinker::deaz
grad_prec * deaz
tinker::pang
real pang
tinker::cang
real cang
tinker::deax
grad_prec * deax
tinker::Angle::IN_PLANE
@ IN_PLANE
tinker::Angle::FOURIER
@ FOURIER
tinker::Angle::HARMONIC
@ HARMONIC
tinker::Angle::LINEAR
@ LINEAR
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::_1radian
constexpr real _1radian
Definition: const.h:17
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
Definition: testrt.h:9
tinker::dk_angle
__device__ void dk_angle(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ deax, grad_prec *__restrict__ deay, grad_prec *__restrict__ deaz, const Angle *__restrict__ angtyp, real angunit, int i, const int(*__restrict__ iang)[4], const real *__restrict__ anat, const real *__restrict__ ak, const real *__restrict__ afld, real cang, real qang, real pang, real sang, const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: angle.h:29
Generated by doxygen 1.9.5