4#include "math/libfunc.h"
25 constexpr bool do_e = Ver::e;
26 constexpr bool do_g = Ver::g;
27 constexpr bool do_v = Ver::v;
31 int ja1 = igrp[ia][0];
32 int ja2 = igrp[ia][1];
33 int jb1 = igrp[ib][0];
34 int jb2 = igrp[ib][1];
39 real weigha = REAL_MAX((
real)1, grpmass[ia]);
41 real xa1 =
x[ka1], ya1 =
y[ka1], za1 =
z[ka1];
42 real xacm = xa1 * weigha, yacm = ya1 * weigha, zacm = za1 * weigha;
44 for (
int j = ja1 + 1; j < ja2; ++j) {
50 if (molec[ka1] != molec[k])
image(xr, yr, zr);
55 weigha = REAL_RECIP(weigha);
57 real weighb = REAL_MAX((
real)1, grpmass[ib]);
59 real xb1 =
x[kb1], yb1 =
y[kb1], zb1 =
z[kb1];
60 real xbcm = xb1 * weighb, ybcm = yb1 * weighb, zbcm = zb1 * weighb;
62 for (
int j = jb1 + 1; j < jb2; ++j) {
68 if (molec[kb1] != molec[k])
image(xr, yr, zr);
73 weighb = REAL_RECIP(weighb);
75 real xr = xacm * weigha - xbcm * weighb;
76 real yr = yacm * weigha - ybcm * weighb;
77 real zr = zacm * weigha - zbcm * weighb;
79 bool intermol = molec[kgrp[ja1]] != molec[kgrp[jb1]];
80 if (intermol)
image(xr, yr, zr);
82 real r = REAL_SQRT(xr * xr + yr * yr + zr * zr);
86 real target = (r < gf1 ? gf1 : (r > gf2 ? gf2 : r));
94 real rinv = (r == 0 ? 1 : REAL_RECIP(r));
95 real de = 2 * force * dt * rinv;
101 for (
int j = ja1; j < ja2; ++j) {
109 for (
int j = jb1; j < jb2; ++j) {
127#pragma acc routine seq
141 constexpr bool do_e = Ver::e;
142 constexpr bool do_g = Ver::g;
143 constexpr bool do_v = Ver::v;
145 int ia =
idfix[i][0];
146 int ib =
idfix[i][1];
154 bool intermol = molec[ia] != molec[ib];
155 if (intermol)
image(xr, yr, zr);
156 real r = REAL_SQRT(xr * xr + yr * yr + zr * zr);
158 if (r < df1) target = df1;
159 if (r > df2) target = df2;
160 real dt = r - target;
167 real rinv = (r == 0 ? 1 : REAL_RECIP(r));
168 real de = 2 * force * dt * rinv;
189#pragma acc routine seq
202 constexpr bool do_e = Ver::e;
203 constexpr bool do_g = Ver::g;
204 constexpr bool do_v = Ver::v;
206 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
208 int ia =
iafix[i][0];
209 int ib =
iafix[i][1];
210 int ic =
iafix[i][2];
225 real xab = xia - xib;
226 real yab = yia - yib;
227 real zab = zia - zib;
228 real xcb = xic - xib;
229 real ycb = yic - yib;
230 real zcb = zic - zib;
231 real rab2 = xab * xab + yab * yab + zab * zab;
232 rab2 = REAL_MAX(rab2, (
real)0.0001);
233 real rcb2 = xcb * xcb + ycb * ycb + zcb * zcb;
234 rcb2 = REAL_MAX(rcb2, (
real)0.0001);
236 real xp = ycb * zab - zcb * yab;
237 real yp = zcb * xab - xcb * zab;
238 real zp = xcb * yab - ycb * xab;
239 real rp = REAL_SQRT(xp * xp + yp * yp + zp * zp);
240 rp = REAL_MAX(rp, (
real)0.0001);
241 real dot = xab * xcb + yab * ycb + zab * zcb;
242 real cosine = dot * REAL_RSQRT(rab2 * rcb2);
243 cosine = REAL_MIN((
real)1, REAL_MAX((
real)-1, cosine));
247 if (angle < af1) target = af1;
248 if (angle > af2) target = af2;
256 real deddt = 2 * force * dt;
257 real terma = -deddt * REAL_RECIP(rab2 * rp);
258 real termc = deddt * REAL_RECIP(rcb2 * rp);
259 real dedxia = terma * (yab * zp - zab * yp);
260 real dedyia = terma * (zab * xp - xab * zp);
261 real dedzia = terma * (xab * yp - yab * xp);
262 real dedxic = termc * (ycb * zp - zcb * yp);
263 real dedyic = termc * (zcb * xp - xcb * zp);
264 real dedzic = termc * (xcb * yp - ycb * xp);
265 real dedxib = -dedxia - dedxic;
266 real dedyib = -dedyia - dedyic;
267 real dedzib = -dedzia - dedzic;
279 vxx = xab * dedxia + xcb * dedxic;
280 vyx = yab * dedxia + ycb * dedxic;
281 vzx = zab * dedxia + zcb * dedxic;
282 vyy = yab * dedyia + ycb * dedyic;
283 vzy = zab * dedyia + zcb * dedyic;
284 vzz = zab * dedzia + zcb * dedzic;
289#pragma acc routine seq
302 constexpr bool do_e = Ver::e;
303 constexpr bool do_g = Ver::g;
304 constexpr bool do_v = Ver::v;
306 if CONSTEXPR (do_v) vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
308 int ia =
itfix[i][0];
309 int ib =
itfix[i][1];
310 int ic =
itfix[i][2];
311 int id =
itfix[i][3];
328 real xba = xib - xia;
329 real yba = yib - yia;
330 real zba = zib - zia;
331 real xcb = xic - xib;
332 real ycb = yic - yib;
333 real zcb = zic - zib;
334 real rcb = REAL_SQRT(xcb * xcb + ycb * ycb + zcb * zcb);
335 rcb = REAL_MAX(rcb, (
real)0.0001);
336 real xdc = xid - xic;
337 real ydc = yid - yic;
338 real zdc = zid - zic;
340 real xt = yba * zcb - ycb * zba;
341 real yt = zba * xcb - zcb * xba;
342 real zt = xba * ycb - xcb * yba;
343 real xu = ycb * zdc - ydc * zcb;
344 real yu = zcb * xdc - zdc * xcb;
345 real zu = xcb * ydc - xdc * ycb;
346 real xtu = yt * zu - yu * zt;
347 real ytu = zt * xu - zu * xt;
348 real ztu = xt * yu - xu * yt;
349 real rt2 = xt * xt + yt * yt + zt * zt;
350 rt2 = REAL_MAX(rt2, (
real)0.0001);
351 real ru2 = xu * xu + yu * yu + zu * zu;
352 ru2 = REAL_MAX(ru2, (
real)0.0001);
353 real rtru = REAL_SQRT(rt2 * ru2);
355 real cosine = (xt * xu + yt * yu + zt * zu) * REAL_RECIP(rtru);
356 real sine = (xcb * xtu + ycb * ytu + zcb * ztu) * REAL_RECIP(rcb * rtru);
357 cosine = REAL_MIN((
real)1.0, REAL_MAX((
real)-1.0, cosine));
359 if (sine < 0) angle = -angle;
361 if (angle > tf1 and angle < tf2)
363 else if (angle > tf1 and tf1 > tf2)
365 else if (angle < tf2 and tf1 > tf2)
368 real t1 = angle - tf1;
369 real t2 = angle - tf2;
370 if (t1 > 180) t1 -= 360;
371 if (t1 < -180) t1 += 360;
372 if (t2 > 180) t2 -= 360;
373 if (t2 < -180) t2 += 360;
374 if (REAL_ABS(t1) < REAL_ABS(t2))
379 real dt = angle - target;
380 if (dt > 180) dt -= 360;
381 if (dt < -180) dt += 360;
389 real dedphi = 2 * force * dt;
391 real xca = xic - xia;
392 real yca = yic - yia;
393 real zca = zic - zia;
394 real xdb = xid - xib;
395 real ydb = yid - yib;
396 real zdb = zid - zib;
398 real rt_inv = REAL_RECIP(rt2 * rcb);
399 real ru_inv = REAL_RECIP(ru2 * rcb);
400 real dedxt = dedphi * (yt * zcb - ycb * zt) * rt_inv;
401 real dedyt = dedphi * (zt * xcb - zcb * xt) * rt_inv;
402 real dedzt = dedphi * (xt * ycb - xcb * yt) * rt_inv;
403 real dedxu = -dedphi * (yu * zcb - ycb * zu) * ru_inv;
404 real dedyu = -dedphi * (zu * xcb - zcb * xu) * ru_inv;
405 real dedzu = -dedphi * (xu * ycb - xcb * yu) * ru_inv;
407 real dedxia = zcb * dedyt - ycb * dedzt;
408 real dedyia = xcb * dedzt - zcb * dedxt;
409 real dedzia = ycb * dedxt - xcb * dedyt;
410 real dedxib = yca * dedzt - zca * dedyt + zdc * dedyu - ydc * dedzu;
411 real dedyib = zca * dedxt - xca * dedzt + xdc * dedzu - zdc * dedxu;
412 real dedzib = xca * dedyt - yca * dedxt + ydc * dedxu - xdc * dedyu;
413 real dedxic = zba * dedyt - yba * dedzt + ydb * dedzu - zdb * dedyu;
414 real dedyic = xba * dedzt - zba * dedxt + zdb * dedxu - xdb * dedzu;
415 real dedzic = yba * dedxt - xba * dedyt + xdb * dedyu - ydb * dedxu;
416 real dedxid = zcb * dedyu - ycb * dedzu;
417 real dedyid = xcb * dedzu - zcb * dedxu;
418 real dedzid = ycb * dedxu - xcb * dedyu;
434 vxx = xcb * (dedxic + dedxid) - xba * dedxia + xdc * dedxid;
435 vyx = ycb * (dedxic + dedxid) - yba * dedxia + ydc * dedxid;
436 vzx = zcb * (dedxic + dedxid) - zba * dedxia + zdc * dedxid;
437 vyy = ycb * (dedyic + dedyid) - yba * dedyia + ydc * dedyid;
438 vzy = zcb * (dedyic + dedyid) - zba * dedyia + zdc * dedyid;
439 vzz = zcb * (dedzic + dedzid) - zba * dedzia + zdc * dedzid;
444#pragma acc routine seq
460 constexpr bool do_e = Ver::e;
461 constexpr bool do_g = Ver::g;
462 constexpr bool do_v = Ver::v;
467 real xr = 0, yr = 0, zr = 0;
472 real r = REAL_SQRT(xr * xr + yr * yr + zr * zr);
473 real dt = REAL_MAX((
real)0, r - radius);
480 real rinv = (r == 0 ? 1 : REAL_RECIP(r));
481 real de = 2 * force * dt * rinv;
void atomic_add(T value, T *buffer, size_t offset=0)
Definition: acc/adddef.h:14
#define SEQ_CUDA
Definition: acc/seqdef.h:12
#define image(x, y, z)
Definition: image.h:210
#define TINKER_IMAGE_PARAMS
Definition: box.h:110
#define restrict
Definition: macro.h:51
#define CONSTEXPR
Definition: macro.h:61
double * mass
Atomic mass.
real * z
Current coordinates used in energy evaluation and neighbor lists.
real * x
Number of the trajectory frames.
real * y
Current coordinates used in energy evaluation and neighbor lists.
int(* igfix)[2]
Group numbers defining each group distance restraint.
real(* gfix)[3]
Force constant and target range for each group distance.
constexpr real _1radian
Definition: const.h:17
constexpr real radian
Definition: const.h:14
float real
Definition: precision.h:80
fixed grad_prec
Definition: precision.h:103
__device__ void dk_geom_group(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ degx, grad_prec *__restrict__ degy, grad_prec *__restrict__ degz, int i, const int(*__restrict__ igfix)[2], const real(*__restrict__ gfix)[3], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z, const double *__restrict__ mass, const int *__restrict__ molec, const int(*__restrict__ igrp)[2], const int *__restrict__ kgrp, const double *__restrict__ grpmass, BoxShape box_shape, real3 lvec1, real3 lvec2, real3 lvec3, real3 recipa, real3 recipb, real3 recipc)
Definition: geom.h:12
__device__ void dk_geom_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__ degx, grad_prec *__restrict__ degy, grad_prec *__restrict__ degz, int i, const int(*__restrict__ iafix)[3], const real(*__restrict__ afix)[3], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: geom.h:192
__device__ void dk_geom_distance(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ degx, grad_prec *__restrict__ degy, grad_prec *__restrict__ degz, int i, const int(*__restrict__ idfix)[2], const real(*__restrict__ dfix)[3], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z, const int *__restrict__ molec, BoxShape box_shape, real3 lvec1, real3 lvec2, real3 lvec3, real3 recipa, real3 recipb, real3 recipc)
Definition: geom.h:130
__device__ void dk_geom_position(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ degx, grad_prec *__restrict__ degy, grad_prec *__restrict__ degz, int i, const int *__restrict__ ipfix, const int(*__restrict__ kpfix)[3], const real *__restrict__ xpfix, const real *__restrict__ ypfix, const real *__restrict__ zpfix, const real(*__restrict__ pfix)[2], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z, BoxShape box_shape, real3 lvec1, real3 lvec2, real3 lvec3, real3 recipa, real3 recipb, real3 recipc)
Definition: geom.h:447
__device__ void dk_geom_torsion(real &__restrict__ e, real &__restrict__ vxx, real &__restrict__ vyx, real &__restrict__ vzx, real &__restrict__ vyy, real &__restrict__ vzy, real &__restrict__ vzz, grad_prec *__restrict__ degx, grad_prec *__restrict__ degy, grad_prec *__restrict__ degz, int i, const int(*__restrict__ itfix)[4], const real(*__restrict__ tfix)[3], const real *__restrict__ x, const real *__restrict__ y, const real *__restrict__ z)
Definition: geom.h:292