Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
opbend.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 {
22#pragma acc routine seq
23template <class Ver>
24SEQ_CUDA
25void dk_opbend(real& restrict e, real& restrict vxx, real& restrict vyx,
26 real& restrict vzx, real& restrict vyy, real& restrict vzy,
27 real& restrict vzz,
28
29 grad_prec* restrict deopbx, grad_prec* restrict deopby,
30 grad_prec* restrict deopbz,
31
32 OPBend opbtyp, real opbunit, int iopbend, const int* restrict iopb,
33 const real* restrict opbk, const int (*restrict iang)[4], real copb,
34 real qopb, real popb, real sopb,
35
36 const real* restrict x, const real* restrict y, const real* restrict z)
37{
38 constexpr bool do_e = Ver::e;
39 constexpr bool do_g = Ver::g;
40 constexpr bool do_v = Ver::v;
41 if CONSTEXPR (do_e) e = 0;
42 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
43
44 const real force = opbk[iopbend];
45 const int i = iopb[iopbend];
46 const int ia = iang[i][0];
47 const int ib = iang[i][1];
48 const int ic = iang[i][2];
49 const int id = iang[i][3];
50
51 real xia = x[ia];
52 real yia = y[ia];
53 real zia = z[ia];
54 real xib = x[ib];
55 real yib = y[ib];
56 real zib = z[ib];
57 real xic = x[ic];
58 real yic = y[ic];
59 real zic = z[ic];
60 real xid = x[id];
61 real yid = y[id];
62 real zid = z[id];
63
64 real xab = xia - xib;
65 real yab = yia - yib;
66 real zab = zia - zib;
67 real xcb = xic - xib;
68 real ycb = yic - yib;
69 real zcb = zic - zib;
70 real xdb = xid - xib;
71 real ydb = yid - yib;
72 real zdb = zid - zib;
73 real xad = xia - xid;
74 real yad = yia - yid;
75 real zad = zia - zid;
76 real xcd = xic - xid;
77 real ycd = yic - yid;
78 real zcd = zic - zid;
79
80 real rab2, rad2;
81 real rcb2, rcd2;
82 real dot;
83 real cc;
84 if (opbtyp == OPBend::WDC) {
85
86 // W-D-C angle between A-B-C plane and B-D vector for D-B<AC
87
88 rab2 = xab * xab + yab * yab + zab * zab;
89 rcb2 = xcb * xcb + ycb * ycb + zcb * zcb;
90 cc = rab2 * rcb2
91 - (xab * xcb + yab * ycb + zab * zcb)
92 * (xab * xcb + yab * ycb + zab * zcb);
93 if CONSTEXPR (do_g) dot = xab * xcb + yab * ycb + zab * zcb;
94 } else if (opbtyp == OPBend::ALLINGER) {
95
96 // Allinger angle between A-C-D plane and D-B vector for D-B<AC
97
98 rad2 = xad * xad + yad * yad + zad * zad;
99 rcd2 = xcd * xcd + ycd * ycd + zcd * zcd;
100 cc = rad2 * rcd2
101 - (xad * xcd + yad * ycd + zad * zcd)
102 * (xad * xcd + yad * ycd + zad * zcd);
103 if CONSTEXPR (do_g) dot = xad * xcd + yad * ycd + zad * zcd;
104 }
105
106 // find the out-of-plane angle bending energy
107
108 real ee = xdb * (yab * zcb - zab * ycb) + ydb * (zab * xcb - xab * zcb)
109 + zdb * (xab * ycb - yab * xcb);
110 real rdb2 = xdb * xdb + ydb * ydb + zdb * zdb;
111 rdb2 = REAL_MAX(rdb2, (real)0.0001);
112
113 if (cc != 0) {
114 real sine = REAL_ABS(ee) * REAL_RSQRT(cc * rdb2);
115 sine = REAL_MIN((real)1, sine);
116 real angle = radian * REAL_ASIN(sine);
117 real dt = angle;
118 real dt2 = dt * dt;
119 real dt3 = dt2 * dt;
120 real dt4 = dt2 * dt2;
121
122 if CONSTEXPR (do_e)
123 e = opbunit * force * dt2
124 * (1 + copb * dt + qopb * dt2 + popb * dt3 + sopb * dt4);
125
126 if CONSTEXPR (do_g) {
127 real deddt = opbunit * force * dt * radian
128 * (2 + 3 * copb * dt + 4 * qopb * dt2 + 5 * popb * dt3
129 + 6 * sopb * dt4);
130 real dedcos = -deddt * REAL_SIGN((real)1, ee)
131 * REAL_RSQRT(cc * rdb2 - ee * ee);
132 real term = ee * REAL_RECIP(cc);
133 real dccdxia, dccdyia, dccdzia;
134 real dccdxic, dccdyic, dccdzic;
135 real dccdxid, dccdyid, dccdzid;
136 if (opbtyp == OPBend::WDC) {
137 dccdxia = (xab * rcb2 - xcb * dot) * term;
138 dccdyia = (yab * rcb2 - ycb * dot) * term;
139 dccdzia = (zab * rcb2 - zcb * dot) * term;
140 dccdxic = (xcb * rab2 - xab * dot) * term;
141 dccdyic = (ycb * rab2 - yab * dot) * term;
142 dccdzic = (zcb * rab2 - zab * dot) * term;
143 dccdxid = 0;
144 dccdyid = 0;
145 dccdzid = 0;
146 } else if (opbtyp == OPBend::ALLINGER) {
147 dccdxia = (xad * rcd2 - xcd * dot) * term;
148 dccdyia = (yad * rcd2 - ycd * dot) * term;
149 dccdzia = (zad * rcd2 - zcd * dot) * term;
150 dccdxic = (xcd * rad2 - xad * dot) * term;
151 dccdyic = (ycd * rad2 - yad * dot) * term;
152 dccdzic = (zcd * rad2 - zad * dot) * term;
153 dccdxid = -dccdxia - dccdxic;
154 dccdyid = -dccdyia - dccdyic;
155 dccdzid = -dccdzia - dccdzic;
156 }
157
158 term = ee * REAL_RECIP(rdb2);
159 real deedxia = ydb * zcb - zdb * ycb;
160 real deedyia = zdb * xcb - xdb * zcb;
161 real deedzia = xdb * ycb - ydb * xcb;
162 real deedxic = yab * zdb - zab * ydb;
163 real deedyic = zab * xdb - xab * zdb;
164 real deedzic = xab * ydb - yab * xdb;
165 real deedxid = ycb * zab - zcb * yab + xdb * term;
166 real deedyid = zcb * xab - xcb * zab + ydb * term;
167 real deedzid = xcb * yab - ycb * xab + zdb * term;
168
169 // compute first derivative components for this angle
170
171 real dedxia = dedcos * (dccdxia + deedxia);
172 real dedyia = dedcos * (dccdyia + deedyia);
173 real dedzia = dedcos * (dccdzia + deedzia);
174 real dedxic = dedcos * (dccdxic + deedxic);
175 real dedyic = dedcos * (dccdyic + deedyic);
176 real dedzic = dedcos * (dccdzic + deedzic);
177 real dedxid = dedcos * (dccdxid + deedxid);
178 real dedyid = dedcos * (dccdyid + deedyid);
179 real dedzid = dedcos * (dccdzid + deedzid);
180 real dedxib = -dedxia - dedxic - dedxid;
181 real dedyib = -dedyia - dedyic - dedyid;
182 real dedzib = -dedzia - dedzic - dedzid;
183
184 atomic_add(dedxia, deopbx, ia);
185 atomic_add(dedyia, deopby, ia);
186 atomic_add(dedzia, deopbz, ia);
187 atomic_add(dedxib, deopbx, ib);
188 atomic_add(dedyib, deopby, ib);
189 atomic_add(dedzib, deopbz, ib);
190 atomic_add(dedxic, deopbx, ic);
191 atomic_add(dedyic, deopby, ic);
192 atomic_add(dedzic, deopbz, ic);
193 atomic_add(dedxid, deopbx, id);
194 atomic_add(dedyid, deopby, id);
195 atomic_add(dedzid, deopbz, id);
196
197 if CONSTEXPR (do_v) {
198 vxx = xab * dedxia + xcb * dedxic + xdb * dedxid;
199 vyx = yab * dedxia + ycb * dedxic + ydb * dedxid;
200 vzx = zab * dedxia + zcb * dedxic + zdb * dedxid;
201 vyy = yab * dedyia + ycb * dedyic + ydb * dedyid;
202 vzy = zab * dedyia + zcb * dedyic + zdb * dedyid;
203 vzz = zab * dedzia + zcb * dedzic + zdb * dedzid;
204 }
205 }
206 }
207}
208}
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::deopby
grad_prec * deopby
tinker::opbk
real * opbk
tinker::OPBend
OPBend
Definition: evalence.h:76
tinker::sopb
real sopb
tinker::popb
real popb
tinker::deopbz
grad_prec * deopbz
tinker::iang
int(* iang)[4]
tinker::copb
real copb
tinker::opbunit
real opbunit
tinker::opbtyp
OPBend opbtyp
tinker::deopbx
grad_prec * deopbx
tinker::iopb
int * iopb
tinker::qopb
real qopb
tinker::OPBend::ALLINGER
@ ALLINGER
tinker::OPBend::WDC
@ WDC
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
Definition: testrt.h:9
tinker::dk_opbend
__device__ void dk_opbend(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ deopbx, grad_prec *__restrict__ deopby, grad_prec *__restrict__ deopbz, OPBend opbtyp, real opbunit, int iopbend, const int *__restrict__ iopb, const real *__restrict__ opbk, const int(*__restrict__ iang)[4], real copb, real qopb, real popb, real sopb, const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: opbend.h:25
Generated by doxygen 1.9.5