Tinker9 70bd052 (Thu Nov 9 12:11:35 2023 -0800)
Loading...
Searching...
No Matches
settle.h
1#pragma once
2#include "ff/precision.h"
3#include "seq/seq.h"
4#include <tinker/detail/units.hh>
5
6namespace tinker {
7#pragma acc routine seq
8template <class HTYPE>
9SEQ_CUDA
10void dk_settle1(time_prec dt, int iw, const pos_prec* restrict xold,
11 const pos_prec* restrict yold, const pos_prec* restrict zold,
12 pos_prec* restrict xnew, pos_prec* restrict ynew, pos_prec* restrict znew,
13 vel_prec* restrict vx, vel_prec* restrict vy, vel_prec* restrict vz,
14 const double* restrict mass, const int (*restrict iratwt)[3],
15 const pos_prec (*restrict kratwt)[3])
16{
17 // atoms a, b, c; lengths ab, ac, bc
18 int ia, ib, ic;
19 double lab, lac, lbc;
20 double m0, m1, m2, invm;
21 ia = iratwt[iw][0];
22 ib = iratwt[iw][1];
23 ic = iratwt[iw][2];
24 lab = kratwt[iw][0];
25 lac = kratwt[iw][1];
26 lbc = kratwt[iw][2];
27 m0 = mass[ia];
28 m1 = mass[ib];
29 m2 = mass[ic];
30 invm = 1 / (m0 + m1 + m2);
31
32 // ABC0 is the triangle at t0.
33 // ABC1 is the triangle at t0+dt in the absence of constraints.
34
35 // global frame vectors AB0 and AC0
36 double xb0, yb0, zb0, xc0, yc0, zc0;
37 xb0 = xold[ib] - xold[ia];
38 yb0 = yold[ib] - yold[ia];
39 zb0 = zold[ib] - zold[ia];
40 xc0 = xold[ic] - xold[ia];
41 yc0 = yold[ic] - yold[ia];
42 zc0 = zold[ic] - zold[ia];
43
44 // global frame centroid D1
45 double xcom, ycom, zcom;
46 xcom = (m0 * xnew[ia] + m1 * xnew[ib] + m2 * xnew[ic]) * invm;
47 ycom = (m0 * ynew[ia] + m1 * ynew[ib] + m2 * ynew[ic]) * invm;
48 zcom = (m0 * znew[ia] + m1 * znew[ib] + m2 * znew[ic]) * invm;
49
50 // global frame vectors DA1, DB1, DC1
51 double xa1, ya1, za1, xb1, yb1, zb1, xc1, yc1, zc1;
52 xa1 = xnew[ia] - xcom;
53 ya1 = ynew[ia] - ycom;
54 za1 = znew[ia] - zcom;
55 xb1 = xnew[ib] - xcom;
56 yb1 = ynew[ib] - ycom;
57 zb1 = znew[ib] - zcom;
58 xc1 = xnew[ic] - xcom;
59 yc1 = ynew[ic] - ycom;
60 zc1 = znew[ic] - zcom;
61
62 // local frame unit vectors x', y', z' and rotation matrix
63 double xakszd, yakszd, zakszd;
64 double xaksxd, yaksxd, zaksxd;
65 double xaksyd, yaksyd, zaksyd;
66 double axlng, aylng, azlng;
67 double trns11, trns21, trns31;
68 double trns12, trns22, trns32;
69 double trns13, trns23, trns33;
70 // z' = AB0 cross AC0
71 xakszd = yb0 * zc0 - zb0 * yc0;
72 yakszd = zb0 * xc0 - xb0 * zc0;
73 zakszd = xb0 * yc0 - yb0 * xc0;
74 // x' = DA1 cross z'
75 xaksxd = ya1 * zakszd - za1 * yakszd;
76 yaksxd = za1 * xakszd - xa1 * zakszd;
77 zaksxd = xa1 * yakszd - ya1 * xakszd;
78 // y' = z' cross x'
79 xaksyd = yakszd * zaksxd - zakszd * yaksxd;
80 yaksyd = zakszd * xaksxd - xakszd * zaksxd;
81 zaksyd = xakszd * yaksxd - yakszd * xaksxd;
82 // normalize x', y', z' to get rotation matrix
83 axlng = 1 / sqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
84 aylng = 1 / sqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
85 azlng = 1 / sqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
86 trns11 = xaksxd * axlng;
87 trns21 = yaksxd * axlng;
88 trns31 = zaksxd * axlng;
89 trns12 = xaksyd * aylng;
90 trns22 = yaksyd * aylng;
91 trns32 = zaksyd * aylng;
92 trns13 = xakszd * azlng;
93 trns23 = yakszd * azlng;
94 trns33 = zakszd * azlng;
95
96 // local frame vectors AB0 and AC0
97 double xb0d, yb0d;
98 double xc0d, yc0d;
99 xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
100 yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
101 xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
102 yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
103
104 // local frame vectors DA1, DB1, DC1
105 // double xa1d, ya1d;
106 double za1d;
107 double xb1d, yb1d, zb1d;
108 double xc1d, yc1d, zc1d;
109 // xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
110 // ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
111 za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
112 xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
113 yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
114 zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
115 xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
116 yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
117 zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
118
119 //====================================================================//
120 // analogous to Eq.A1 for arbitrary triangular constraint
121 // local frame coordinates of abc0'
122 // a0' = (0,yra,0); b0' = (xrb,yrb,0); c0' = (xrc,yrc,0)
123 // In Eq.A1, ra, rb, and rc are lengths.
124 // Here we use the SIGNED coordinates.
125 double yra, xrb, yrb, xrc, yrc;
126 {
127 // first let A = (0,0), B = (lab,0), C = lac(cosA,sinA)
128 // centroid G = (xg,yg)
129 double cosA, sinA, xg, yg;
130 cosA = (lab * lab + lac * lac - lbc * lbc) / (2 * lab * lac);
131 sinA = sqrt(1 - cosA * cosA);
132 xg = (m1 * lab + m2 * lac * cosA) * invm;
133 yg = m2 * lac * sinA * invm;
134
135 // vectors GA, GB, GC
136 double xga, yga, xgb, ygb, xgc, ygc;
137 double lga, lgb, lgc;
138 double cosAGB, cosAGC, sinAGB, sinAGC;
139 xga = -xg;
140 yga = -yg;
141 xgb = lab - xg;
142 ygb = -yg;
143 xgc = lac * cosA - xg;
144 ygc = lac * sinA - yg;
145 lga = sqrt(xga * xga + yga * yga);
146 lgb = sqrt(xgb * xgb + ygb * ygb);
147 lgc = sqrt(xgc * xgc + ygc * ygc);
148 cosAGB = (xga * xgb + yga * ygb) / (lga * lgb);
149 cosAGC = (xga * xgc + yga * ygc) / (lga * lgc);
150 sinAGB = sqrt(1 - cosAGB * cosAGB);
151 sinAGC = sqrt(1 - cosAGC * cosAGC);
152
153 yra = lga;
154 xrb = -sinAGB * lgb;
155 yrb = cosAGB * lgb;
156 xrc = sinAGC * lgc;
157 yrc = cosAGC * lgc;
158 }
159
160 //====================================================================//
161 // Eq.A5
162 double sinphi, cosphi, sinpsi, cospsi;
163 sinphi = za1d / yra;
164 cosphi = sqrt(1 - sinphi * sinphi);
165 sinpsi = ((zb1d - zc1d) - (yrb - yrc) * sinphi) / ((xrc - xrb) * cosphi);
166 cospsi = sqrt(1 - sinpsi * sinpsi);
167
168 //====================================================================//
169 // Eq.A3 local frame vectors da2', db2', dc2'
170 double ya2d;
171 double xb2d, yb2d;
172 double xc2d, yc2d;
173 ya2d = yra * cosphi;
174 xb2d = xrb * cospsi;
175 yb2d = yrb * cosphi + xrb * sinphi * sinpsi;
176 xc2d = xrc * cospsi;
177 yc2d = yrc * cosphi + xrc * sinphi * sinpsi;
178 // adjust numerical error for xb2' and xc2'
179 // deltx**2 + ybc**2 + zbc**2 = lbc**2, where ybc**2 + zbc**2 = hh2
180 double deltx, hh2;
181 hh2 = (yc2d - yb2d) * (yc2d - yb2d) + (zc1d - zb1d) * (zc1d - zb1d);
182 // adjusted xbc
183 deltx = sqrt(lbc * lbc - hh2);
184 // adjustment deltx - (xc - xb)
185 deltx = deltx - xc2d + xb2d;
186 // m0 xa(=0) + m1 xb + m2 xc = 0
187 xb2d -= deltx * m2 / (m1 + m2);
188 xc2d += deltx * m1 / (m1 + m2);
189
190 //====================================================================//
191 // Eq.A15
192 double alpa, beta, gama, al2be2, sinthe, costhe;
193 alpa = m1 * (yb2d * yb0d + xb2d * xb0d);
194 alpa += m2 * (yc2d * yc0d + xc2d * xc0d);
195 beta = m1 * (yb2d * xb0d - xb2d * yb0d);
196 beta += m2 * (yc2d * xc0d - xc2d * yc0d);
197 gama = m1 * (yb1d * xb0d - xb1d * yb0d);
198 gama += m2 * (yc1d * xc0d - xc1d * yc0d);
199 al2be2 = alpa * alpa + beta * beta;
200 sinthe = (alpa * gama - beta * sqrt(al2be2 - gama * gama)) / al2be2;
201 costhe = sqrt(1 - sinthe * sinthe);
202
203 //====================================================================//
204 // Eq.A4 local frame vectors da3', db3', dc3'
205 double xa3d, ya3d, za3d;
206 double xb3d, yb3d, zb3d;
207 double xc3d, yc3d, zc3d;
208 xa3d = -ya2d * sinthe;
209 ya3d = ya2d * costhe;
210 za3d = za1d;
211 xb3d = xb2d * costhe - yb2d * sinthe;
212 yb3d = xb2d * sinthe + yb2d * costhe;
213 zb3d = zb1d;
214 xc3d = xc2d * costhe - yc2d * sinthe;
215 yc3d = xc2d * sinthe + yc2d * costhe;
216 zc3d = zc1d;
217
218 // global frame coordinates DA3, DB3, DC3
219 double xa3, ya3, za3;
220 double xb3, yb3, zb3;
221 double xc3, yc3, zc3;
222 xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
223 ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
224 za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
225 xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
226 yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
227 zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
228 xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
229 yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
230 zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
231
232 xnew[ia] = xcom + xa3;
233 ynew[ia] = ycom + ya3;
234 znew[ia] = zcom + za3;
235 xnew[ib] = xcom + xb3;
236 ynew[ib] = ycom + yb3;
237 znew[ib] = zcom + zb3;
238 xnew[ic] = xcom + xc3;
239 ynew[ic] = ycom + yc3;
240 znew[ic] = zcom + zc3;
241
242 //====================================================================//
243 // velocity correction due to the position constraints
244
245 // This code is not in the original SETTLE paper, but is necessary for
246 // the velocity verlet integrator.
247 if CONSTEXPR (not eq<HTYPE, SHAKE>()) {
248 double invdt = 1 / dt;
249 vx[ia] += (xa3 - xa1) * invdt;
250 vy[ia] += (ya3 - ya1) * invdt;
251 vz[ia] += (za3 - za1) * invdt;
252 vx[ib] += (xb3 - xb1) * invdt;
253 vy[ib] += (yb3 - yb1) * invdt;
254 vz[ib] += (zb3 - zb1) * invdt;
255 vx[ic] += (xc3 - xc1) * invdt;
256 vy[ic] += (yc3 - yc1) * invdt;
257 vz[ic] += (zc3 - zc1) * invdt;
258 }
259}
260
261#pragma acc routine seq
262template <bool DO_V>
263SEQ_CUDA
264void dk_settle2(time_prec dt, int iw, vel_prec* restrict vx,
265 vel_prec* restrict vy, vel_prec* restrict vz, const pos_prec* restrict xpos,
266 const pos_prec* restrict ypos, const pos_prec* restrict zpos,
267 const double* restrict mass, const int (*restrict iratwt)[3],
268 double& restrict vxx, double& restrict vyx, double& restrict vzx,
269 double& restrict vyy, double& restrict vzy, double& restrict vzz)
270{
271 int ia, ib, ic;
272 double m0, m1, m2;
273 ia = iratwt[iw][0];
274 ib = iratwt[iw][1];
275 ic = iratwt[iw][2];
276 m0 = mass[ia];
277 m1 = mass[ib];
278 m2 = mass[ic];
279
280 // vectors AB, BC, CA
281 double xab, yab, zab;
282 double xbc, ybc, zbc;
283 double xca, yca, zca;
284 xab = xpos[ib] - xpos[ia];
285 yab = ypos[ib] - ypos[ia];
286 zab = zpos[ib] - zpos[ia];
287 xbc = xpos[ic] - xpos[ib];
288 ybc = ypos[ic] - ypos[ib];
289 zbc = zpos[ic] - zpos[ib];
290 xca = xpos[ia] - xpos[ic];
291 yca = ypos[ia] - ypos[ic];
292 zca = zpos[ia] - zpos[ic];
293
294 // unit vectors eAB, eBC, eCA
295 double ablng, bclng, calng;
296 double xeab, yeab, zeab;
297 double xebc, yebc, zebc;
298 double xeca, yeca, zeca;
299 ablng = 1 / sqrt(xab * xab + yab * yab + zab * zab);
300 bclng = 1 / sqrt(xbc * xbc + ybc * ybc + zbc * zbc);
301 calng = 1 / sqrt(xca * xca + yca * yca + zca * zca);
302 xeab = xab * ablng;
303 yeab = yab * ablng;
304 zeab = zab * ablng;
305 xebc = xbc * bclng;
306 yebc = ybc * bclng;
307 zebc = zbc * bclng;
308 xeca = xca * calng;
309 yeca = yca * calng;
310 zeca = zca * calng;
311
312 // velocity vectors vAB, vBC, vCA
313 double xvab, yvab, zvab;
314 double xvbc, yvbc, zvbc;
315 double xvca, yvca, zvca;
316 xvab = vx[ib] - vx[ia];
317 yvab = vy[ib] - vy[ia];
318 zvab = vz[ib] - vz[ia];
319 xvbc = vx[ic] - vx[ib];
320 yvbc = vy[ic] - vy[ib];
321 zvbc = vz[ic] - vz[ib];
322 xvca = vx[ia] - vx[ic];
323 yvca = vy[ia] - vy[ic];
324 zvca = vz[ia] - vz[ic];
325
326 double vabab, vbcbc, vcaca;
327 vabab = xvab * xeab + yvab * yeab + zvab * zeab;
328 vbcbc = xvbc * xebc + yvbc * yebc + zvbc * zebc;
329 vcaca = xvca * xeca + yvca * yeca + zvca * zeca;
330
331 double cosa, cosb, cosc;
332 cosa = -xeab * xeca - yeab * yeca - zeab * zeca;
333 cosb = -xebc * xeab - yebc * yeab - zebc * zeab;
334 cosc = -xeca * xebc - yeca * yebc - zeca * zebc;
335
336 //====================================================================//
337 // Eq.B1
338 // M (tab,tbc,tca) = m2v, where
339 // m2v = (ma mb vab0, mb mc vbc0, mc ma vca0),
340 // (a1,a2,a3)
341 // M = (b1,b2,b3)
342 // (c1,c2,c3)
343 double a1, a2, a3, b1, b2, b3, c1, c2, c3, denom;
344 a1 = m0 + m1;
345 a2 = m0 * cosb;
346 a3 = m1 * cosa;
347 b1 = m2 * cosb;
348 b2 = m1 + m2;
349 b3 = m1 * cosc;
350 c1 = m2 * cosa;
351 c2 = m0 * cosc;
352 c3 = m2 + m0;
353 // det(M)
354 denom = a1 * (b2 * c3 - b3 * c2) + a2 * (b3 * c1 - b1 * c3)
355 + a3 * (b1 * c2 - b2 * c1);
356
357 // inverse(M)*det(M)
358 double av1, av2, av3, bv1, bv2, bv3, cv1, cv2, cv3;
359 av1 = b2 * c3 - b3 * c2;
360 av2 = c2 * a3 - c3 * a2;
361 av3 = a2 * b3 - a3 * b2;
362 bv1 = b3 * c1 - b1 * c3;
363 bv2 = c3 * a1 - c1 * a3;
364 bv3 = a3 * b1 - a1 * b3;
365 cv1 = b1 * c2 - b2 * c1;
366 cv2 = c1 * a2 - c2 * a1;
367 cv3 = a1 * b2 - a2 * b1;
368
369 // t = inverse(M)*m2v
370 double tabd, tbcd, tcad;
371 tabd = av1 * m0 * m1 * vabab + av2 * m1 * m2 * vbcbc + av3 * m2 * m0 * vcaca;
372 tbcd = bv1 * m0 * m1 * vabab + bv2 * m1 * m2 * vbcbc + bv3 * m2 * m0 * vcaca;
373 tcad = cv1 * m0 * m1 * vabab + cv2 * m1 * m2 * vbcbc + cv3 * m2 * m0 * vcaca;
374
375 denom = 1 / denom;
376 vx[ia] += (xeab * tabd - xeca * tcad) / m0 * denom;
377 vz[ia] += (zeab * tabd - zeca * tcad) / m0 * denom;
378 vy[ia] += (yeab * tabd - yeca * tcad) / m0 * denom;
379 vx[ib] += (xebc * tbcd - xeab * tabd) / m1 * denom;
380 vy[ib] += (yebc * tbcd - yeab * tabd) / m1 * denom;
381 vz[ib] += (zebc * tbcd - zeab * tabd) / m1 * denom;
382 vx[ic] += (xeca * tcad - xebc * tbcd) / m2 * denom;
383 vy[ic] += (yeca * tcad - yebc * tbcd) / m2 * denom;
384 vz[ic] += (zeca * tcad - zebc * tbcd) / m2 * denom;
385
386 if CONSTEXPR (DO_V) {
387 const double vterm = 2 / (dt * units::ekcal);
388 double xterm, yterm, zterm;
389 vxx = 0, vyx = 0, vzx = 0, vyy = 0, vzy = 0, vzz = 0;
390 xterm = xeab * tabd * denom * vterm;
391 yterm = yeab * tabd * denom * vterm;
392 zterm = zeab * tabd * denom * vterm;
393 vxx += xab * xterm;
394 vyx += yab * xterm;
395 vzx += zab * xterm;
396 vyy += yab * yterm;
397 vzy += zab * yterm;
398 vzz += zab * zterm;
399 xterm = xebc * tbcd * denom * vterm;
400 yterm = yebc * tbcd * denom * vterm;
401 zterm = zebc * tbcd * denom * vterm;
402 vxx += xbc * xterm;
403 vyx += ybc * xterm;
404 vzx += zbc * xterm;
405 vyy += ybc * yterm;
406 vzy += zbc * yterm;
407 vzz += zbc * zterm;
408 xterm = xeca * tcad * denom * vterm;
409 yterm = yeca * tcad * denom * vterm;
410 zterm = zeca * tcad * denom * vterm;
411 vxx += xca * xterm;
412 vyx += yca * xterm;
413 vzx += zca * xterm;
414 vyy += yca * yterm;
415 vzy += zca * yterm;
416 vzz += zca * zterm;
417 }
418}
419}
SEQ_CUDA
#define SEQ_CUDA
Definition: acc/seqdef.h:12
restrict
#define restrict
Definition: macro.h:51
CONSTEXPR
#define CONSTEXPR
Definition: macro.h:61
tinker::mass
double * mass
Atomic mass.
tinker::ypos
pos_prec * ypos
Definition: atom.h:62
tinker::zpos
pos_prec * zpos
Definition: atom.h:62
tinker::xpos
pos_prec * xpos
tinker::iratwt
int(* iratwt)[3]
tinker::kratwt
pos_prec(* kratwt)[3]
tinker::vy
vel_prec * vy
Velocities.
Definition: pq.h:76
tinker::vx
vel_prec * vx
Velocities.
tinker::vz
vel_prec * vz
Velocities.
Definition: pq.h:76
tinker::pos_prec
mixed pos_prec
Floating-point type for coordinates.
Definition: precision.h:93
tinker::time_prec
mixed time_prec
Floating-point type for time.
Definition: precision.h:90
tinker::vel_prec
mixed vel_prec
Floating-point type for velocities.
Definition: precision.h:92
tinker
Definition: testrt.h:9
tinker::dk_settle2
__device__ void dk_settle2(time_prec dt, int iw, vel_prec *__restrict__ vx, vel_prec *__restrict__ vy, vel_prec *__restrict__ vz, const pos_prec *__restrict__ xpos, const pos_prec *__restrict__ ypos, const pos_prec *__restrict__ zpos, const double *__restrict__ mass, const int(*__restrict__ iratwt)[3], double &__restrict__ vxx, double &__restrict__ vyx, double &__restrict__ vzx, double &__restrict__ vyy, double &__restrict__ vzy, double &__restrict__ vzz)
Definition: settle.h:264
tinker::dk_settle1
__device__ void dk_settle1(time_prec dt, int iw, const pos_prec *__restrict__ xold, const pos_prec *__restrict__ yold, const pos_prec *__restrict__ zold, pos_prec *__restrict__ xnew, pos_prec *__restrict__ ynew, pos_prec *__restrict__ znew, vel_prec *__restrict__ vx, vel_prec *__restrict__ vy, vel_prec *__restrict__ vz, const double *__restrict__ mass, const int(*__restrict__ iratwt)[3], const pos_prec(*__restrict__ kratwt)[3])
Definition: settle.h:10
Generated by doxygen 1.9.5