Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
strbnd.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_strbnd(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 debax, grad_prec* restrict debay,
16 grad_prec* restrict debaz,
17
18 real stbnunit, int istrbnd, const int (*restrict isb)[3],
19 const real (*restrict sbk)[2],
20
21 const real* restrict bl, const int (*restrict iang)[4],
22 const real* restrict anat,
23
24 const real* restrict x, const real* restrict y, const real* restrict z)
25{
26 constexpr bool do_e = Ver::e;
27 constexpr bool do_g = Ver::g;
28 constexpr bool do_v = Ver::v;
29 if CONSTEXPR (do_e) e = 0;
30 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
31
32 int i = isb[istrbnd][0];
33 int j = isb[istrbnd][1];
34 int k = isb[istrbnd][2];
35 int ia = iang[i][0];
36 int ib = iang[i][1];
37 int ic = iang[i][2];
38 real force1 = sbk[istrbnd][0];
39 real force2 = sbk[istrbnd][1];
40
41 real xia = x[ia];
42 real yia = y[ia];
43 real zia = z[ia];
44 real xib = x[ib];
45 real yib = y[ib];
46 real zib = z[ib];
47 real xic = x[ic];
48 real yic = y[ic];
49 real zic = z[ic];
50
51 real xab = xia - xib;
52 real yab = yia - yib;
53 real zab = zia - zib;
54 real xcb = xic - xib;
55 real ycb = yic - yib;
56 real zcb = zic - zib;
57
58 real rab2 = xab * xab + yab * yab + zab * zab;
59 real rcb2 = xcb * xcb + ycb * ycb + zcb * zcb;
60
61 if (REAL_MIN(rab2, rcb2) != 0) {
62 real rab = REAL_SQRT(rab2);
63 real rcb = REAL_SQRT(rcb2);
64 real xp = ycb * zab - zcb * yab;
65 real yp = zcb * xab - xcb * zab;
66 real zp = xcb * yab - ycb * xab;
67 real rp = REAL_SQRT(xp * xp + yp * yp + zp * zp);
68 rp = REAL_MAX(rp, (real)(0.0001));
69 real dot = xab * xcb + yab * ycb + zab * zcb;
70 real cosine = dot * REAL_RECIP(rab * rcb);
71 cosine = REAL_MIN((real)1, REAL_MAX((real)-1, cosine));
72 real angle = radian * REAL_ACOS(cosine);
73
74 real dt = angle - anat[i];
75 real term1 = -radian * REAL_RECIP(rab2 * rp);
76 real term2 = radian * REAL_RECIP(rcb2 * rp);
77 real ddtdxia = term1 * (yab * zp - zab * yp);
78 real ddtdyia = term1 * (zab * xp - xab * zp);
79 real ddtdzia = term1 * (xab * yp - yab * xp);
80 real ddtdxic = term2 * (ycb * zp - zcb * yp);
81 real ddtdyic = term2 * (zcb * xp - xcb * zp);
82 real ddtdzic = term2 * (xcb * yp - ycb * xp);
83
84 real dr1 = rab - bl[j];
85 term1 = REAL_RECIP(rab);
86 real dr2 = rcb - bl[k];
87 term2 = REAL_RECIP(rcb);
88 real ddrdxia = term1 * xab;
89 real ddrdyia = term1 * yab;
90 real ddrdzia = term1 * zab;
91 real ddrdxic = term2 * xcb;
92 real ddrdyic = term2 * ycb;
93 real ddrdzic = term2 * zcb;
94
95 term1 = stbnunit * force1;
96 term2 = stbnunit * force2;
97 real termr = term1 * dr1 + term2 * dr2;
98
99 if CONSTEXPR (do_e) e = termr * dt;
100
101 if CONSTEXPR (do_g) {
102 real term1t = term1 * dt;
103 real term2t = term2 * dt;
104 real dedxia = term1t * ddrdxia + termr * ddtdxia;
105 real dedyia = term1t * ddrdyia + termr * ddtdyia;
106 real dedzia = term1t * ddrdzia + termr * ddtdzia;
107 real dedxic = term2t * ddrdxic + termr * ddtdxic;
108 real dedyic = term2t * ddrdyic + termr * ddtdyic;
109 real dedzic = term2t * ddrdzic + termr * ddtdzic;
110 real dedxib = -dedxia - dedxic;
111 real dedyib = -dedyia - dedyic;
112 real dedzib = -dedzia - dedzic;
113
114 atomic_add(dedxia, debax, ia);
115 atomic_add(dedyia, debay, ia);
116 atomic_add(dedzia, debaz, ia);
117 atomic_add(dedxib, debax, ib);
118 atomic_add(dedyib, debay, ib);
119 atomic_add(dedzib, debaz, ib);
120 atomic_add(dedxic, debax, ic);
121 atomic_add(dedyic, debay, ic);
122 atomic_add(dedzic, debaz, ic);
123
124 if CONSTEXPR (do_v) {
125 vxx = xab * dedxia + xcb * dedxic;
126 vyx = yab * dedxia + ycb * dedxic;
127 vzx = zab * dedxia + zcb * dedxic;
128 vyy = yab * dedyia + ycb * dedyic;
129 vzy = zab * dedyia + zcb * dedyic;
130 vzz = zab * dedzia + zcb * dedzic;
131 }
132 }
133 }
134}
135}
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::iang
int(* iang)[4]
tinker::stbnunit
real stbnunit
tinker::bl
real * bl
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::debay
grad_prec * debay
tinker::isb
int(* isb)[3]
tinker::sbk
real(* sbk)[2]
tinker::debax
grad_prec * debax
tinker::debaz
grad_prec * debaz
tinker
Definition: testrt.h:9
tinker::dk_strbnd
__device__ void dk_strbnd(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ debax, grad_prec *__restrict__ debay, grad_prec *__restrict__ debaz, real stbnunit, int istrbnd, const int(*__restrict__ isb)[3], const real(*__restrict__ sbk)[2], const real *__restrict__ bl, const int(*__restrict__ iang)[4], const real *__restrict__ anat, const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: strbnd.h:11
Generated by doxygen 1.9.5