1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROwing Monsters And Cloning Shrimps
44 #include "gmx_fatal.h"
67 typedef struct gmx_settledata
74 static void init_proj_matrix(settleparam_t *p,
75 real invmO, real invmH, real dOH, real dHH)
82 /* We normalize the inverse masses with imO for the matrix inversion.
83 * so we can keep using masses of almost zero for frozen particles,
84 * without running out of the float range in m_inv.
89 /* Construct the constraint coupling matrix */
90 mat[0][0] = imOn + imHn;
91 mat[0][1] = imOn*(1 - 0.5*dHH*dHH/(dOH*dOH));
92 mat[0][2] = imHn*0.5*dHH/dOH;
93 mat[1][1] = mat[0][0];
94 mat[1][2] = mat[0][2];
95 mat[2][2] = imHn + imHn;
96 mat[1][0] = mat[0][1];
97 mat[2][0] = mat[0][2];
98 mat[2][1] = mat[1][2];
100 m_inv(mat, p->invmat);
102 msmul(p->invmat, 1/p->imO, p->invmat);
108 static void settleparam_init(settleparam_t *p,
109 real mO, real mH, real invmO, real invmH,
121 p->ra = 2.0*mH*sqrt(dOH*dOH - p->rc*p->rc)/wohh;
122 p->rb = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
125 /* For projection: connection matrix inversion */
126 init_proj_matrix(p, invmO, invmH, dOH, dHH);
130 fprintf(debug, "wh =%g, rc = %g, ra = %g\n",
131 p->wh, p->rc, p->ra);
132 fprintf(debug, "rb = %g, irc2 = %g, dHH = %g, dOH = %g\n",
133 p->rb, p->irc2, p->dHH, p->dOH);
137 gmx_settledata_t settle_init(real mO, real mH, real invmO, real invmH,
140 gmx_settledata_t settled;
144 settleparam_init(&settled->massw, mO, mH, invmO, invmH, dOH, dHH);
146 settleparam_init(&settled->mass1, 1.0, 1.0, 1.0, 1.0, dOH, dHH);
152 static void check_cons(FILE *fp, char *title, real x[], int OW1, int HW2, int HW3)
154 rvec dOH1, dOH2, dHH;
157 for (m = 0; (m < DIM); m++)
159 dOH1[m] = x[OW1+m]-x[HW2+m];
160 dOH2[m] = x[OW1+m]-x[HW3+m];
161 dHH[m] = x[HW2+m]-x[HW3+m];
163 fprintf(fp, "%10s, OW1=%3d, HW2=%3d, HW3=%3d, dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
164 title, OW1/DIM, HW2/DIM, HW3/DIM, norm(dOH1), norm(dOH2), norm(dHH));
169 void settle_proj(FILE *fp,
170 gmx_settledata_t settled, int econq,
171 int nsettle, t_iatom iatoms[],
174 rvec *der, rvec *derp,
175 int calcvir_atom_end, tensor vir_r_m_dder,
178 /* Settle for projection out constraint components
179 * of derivatives of the coordinates.
180 * Berk Hess 2008-1-10
184 real imO, imH, dOH, dHH, invdOH, invdHH;
186 int i, m, m2, ow1, hw2, hw3;
187 rvec roh2, roh3, rhh, dc, fc, fcv;
188 rvec derm[3], derpm[3];
189 real invvscale, vscale_nhc, veta;
192 calcvir_atom_end *= DIM;
194 if (econq == econqForce)
204 copy_mat(p->invmat, invmat);
210 veta = vetavar->veta;
211 vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
217 for (i = 0; i < nsettle; i++)
224 for (m = 0; m < DIM; m++)
226 /* in the velocity case, these are the velocities, so we
227 need to modify with the pressure control velocities! */
229 derm[0][m] = vscale_nhc*der[ow1][m] + veta*x[ow1][m];
230 derm[1][m] = vscale_nhc*der[hw2][m] + veta*x[hw2][m];
231 derm[2][m] = vscale_nhc*der[hw3][m] + veta*x[hw3][m];
238 rvec_sub(x[ow1], x[hw2], roh2);
239 rvec_sub(x[ow1], x[hw3], roh3);
240 rvec_sub(x[hw2], x[hw3], rhh);
244 pbc_dx_aiuc(pbc, x[ow1], x[hw2], roh2);
245 pbc_dx_aiuc(pbc, x[ow1], x[hw3], roh3);
246 pbc_dx_aiuc(pbc, x[hw2], x[hw3], rhh);
248 svmul(invdOH, roh2, roh2);
249 svmul(invdOH, roh3, roh3);
250 svmul(invdHH, rhh, rhh);
253 /* Determine the projections of der(modified) on the bonds */
255 for (m = 0; m < DIM; m++)
257 dc[0] += (derm[0][m] - derm[1][m])*roh2[m];
258 dc[1] += (derm[0][m] - derm[2][m])*roh3[m];
259 dc[2] += (derm[1][m] - derm[2][m])*rhh [m];
263 /* Determine the correction for the three bonds */
264 mvmul(invmat, dc, fc);
267 /* divide velocity by vscale_nhc for determining constrained velocities, since they
268 have not yet been multiplied */
269 svmul(1.0/vscale_nhc, fc, fcv);
272 /* Subtract the corrections from derp */
273 for (m = 0; m < DIM; m++)
275 derp[ow1][m] -= imO*( fcv[0]*roh2[m] + fcv[1]*roh3[m]);
276 derp[hw2][m] -= imH*(-fcv[0]*roh2[m] + fcv[2]*rhh [m]);
277 derp[hw3][m] -= imH*(-fcv[1]*roh3[m] - fcv[2]*rhh [m]);
282 if (ow1 < calcvir_atom_end)
284 /* Determining r \dot m der is easy,
285 * since fc contains the mass weighted corrections for der.
288 for (m = 0; m < DIM; m++)
290 for (m2 = 0; m2 < DIM; m2++)
292 vir_r_m_dder[m][m2] +=
293 dOH*roh2[m]*roh2[m2]*fcv[0] +
294 dOH*roh3[m]*roh3[m2]*fcv[1] +
295 dHH*rhh [m]*rhh [m2]*fcv[2];
301 if (calcvir_atom_end > 0)
303 /* Correct r_m_dder, which will be used to calcualate the virial;
304 * we need to use the unscaled multipliers in the virial.
306 msmul(vir_r_m_dder, 1.0/vetavar->vscale, vir_r_m_dder);
311 void csettle(gmx_settledata_t settled,
312 int nsettle, t_iatom iatoms[],
314 real b4[], real after[],
315 real invdt, real *v, int CalcVirAtomEnd,
320 /* ***************************************************************** */
322 /* Subroutine : setlep - reset positions of TIP3P waters ** */
323 /* Author : Shuichi Miyamoto ** */
324 /* Date of last update : Oct. 1, 1992 ** */
326 /* Reference for the SETTLE algorithm ** */
327 /* S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
329 /* ***************************************************************** */
331 /* Initialized data */
333 real wh, ra, rb, rc, irc2;
334 real mOs, mHs, invdts;
336 /* Local variables */
337 real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2;
338 real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22,
339 trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe,
340 cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd,
341 xb0, yb0, zb0, xc0, yc0, zc0, xa1;
342 real ya1, za1, xb1, yb1;
343 real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3,
344 xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d,
345 za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d,
346 xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d;
348 real dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
349 real mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
353 int i, ow1, hw2, hw3;
355 rvec dx, sh_hw2 = {0, 0, 0}, sh_hw3 = {0, 0, 0};
370 mOs = p->mO / vetavar->rvscale;
371 mHs = p->mH / vetavar->rvscale;
372 invdts = invdt / vetavar->rscale;
377 for (i = 0; i < nsettle; ++i)
380 /* --- Step1 A1' --- */
381 ow1 = iatoms[i*4+1] * 3;
382 hw2 = iatoms[i*4+2] * 3;
383 hw3 = iatoms[i*4+3] * 3;
386 xb0 = b4[hw2 + XX] - b4[ow1 + XX];
387 yb0 = b4[hw2 + YY] - b4[ow1 + YY];
388 zb0 = b4[hw2 + ZZ] - b4[ow1 + ZZ];
389 xc0 = b4[hw3 + XX] - b4[ow1 + XX];
390 yc0 = b4[hw3 + YY] - b4[ow1 + YY];
391 zc0 = b4[hw3 + ZZ] - b4[ow1 + ZZ];
394 rvec_sub(after+hw2, after+ow1, doh2);
395 rvec_sub(after+hw3, after+ow1, doh3);
400 pbc_dx_aiuc(pbc, b4+hw2, b4+ow1, dx);
404 pbc_dx_aiuc(pbc, b4+hw3, b4+ow1, dx);
409 /* Tedious way of doing pbc */
410 is = pbc_dx_aiuc(pbc, after+hw2, after+ow1, doh2);
417 sh_hw2[XX] = after[hw2 + XX] - (after[ow1 + XX] + doh2[XX]);
418 sh_hw2[YY] = after[hw2 + YY] - (after[ow1 + YY] + doh2[YY]);
419 sh_hw2[ZZ] = after[hw2 + ZZ] - (after[ow1 + ZZ] + doh2[ZZ]);
420 rvec_dec(after+hw2, sh_hw2);
422 is = pbc_dx_aiuc(pbc, after+hw3, after+ow1, doh3);
429 sh_hw3[XX] = after[hw3 + XX] - (after[ow1 + XX] + doh3[XX]);
430 sh_hw3[YY] = after[hw3 + YY] - (after[ow1 + YY] + doh3[YY]);
431 sh_hw3[ZZ] = after[hw3 + ZZ] - (after[ow1 + ZZ] + doh3[ZZ]);
432 rvec_dec(after+hw3, sh_hw3);
436 /* Not calculating the center of mass using the oxygen position
437 * and the O-H distances, as done below, will make SETTLE
438 * the largest source of energy drift for simulations of water,
439 * as then the oxygen coordinate is multiplied by 0.89 at every step,
440 * which can then transfer a systematic rounding to the oxygen velocity.
442 xa1 = -(doh2[XX] + doh3[XX]) * wh;
443 ya1 = -(doh2[YY] + doh3[YY]) * wh;
444 za1 = -(doh2[ZZ] + doh3[ZZ]) * wh;
446 xcom = after[ow1 + XX] - xa1;
447 ycom = after[ow1 + YY] - ya1;
448 zcom = after[ow1 + ZZ] - za1;
450 xb1 = after[hw2 + XX] - xcom;
451 yb1 = after[hw2 + YY] - ycom;
452 zb1 = after[hw2 + ZZ] - zcom;
453 xc1 = after[hw3 + XX] - xcom;
454 yc1 = after[hw3 + YY] - ycom;
455 zc1 = after[hw3 + ZZ] - zcom;
458 xakszd = yb0 * zc0 - zb0 * yc0;
459 yakszd = zb0 * xc0 - xb0 * zc0;
460 zakszd = xb0 * yc0 - yb0 * xc0;
461 xaksxd = ya1 * zakszd - za1 * yakszd;
462 yaksxd = za1 * xakszd - xa1 * zakszd;
463 zaksxd = xa1 * yakszd - ya1 * xakszd;
464 xaksyd = yakszd * zaksxd - zakszd * yaksxd;
465 yaksyd = zakszd * xaksxd - xakszd * zaksxd;
466 zaksyd = xakszd * yaksxd - yakszd * xaksxd;
469 axlng = gmx_invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
470 aylng = gmx_invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
471 azlng = gmx_invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
473 trns11 = xaksxd * axlng;
474 trns21 = yaksxd * axlng;
475 trns31 = zaksxd * axlng;
476 trns12 = xaksyd * aylng;
477 trns22 = yaksyd * aylng;
478 trns32 = zaksyd * aylng;
479 trns13 = xakszd * azlng;
480 trns23 = yakszd * azlng;
481 trns33 = zakszd * azlng;
484 xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
485 yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
486 xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
487 yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
489 xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
490 ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
492 za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
493 xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
494 yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
495 zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
496 xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
497 yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
498 zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
501 sinphi = za1d * gmx_invsqrt(ra*ra);
502 tmp = 1.0 - sinphi * sinphi;
509 tmp2 = gmx_invsqrt(tmp);
511 sinpsi = (zb1d - zc1d) * irc2 * tmp2;
512 tmp2 = 1.0 - sinpsi * sinpsi;
519 cospsi = tmp2*gmx_invsqrt(tmp2);
529 t2 = rc * sinpsi * sinphi;
534 /* --- Step3 al,be,ga --- */
535 alpa = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
536 beta = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
537 gama = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
538 al2be2 = alpa * alpa + beta * beta;
539 tmp2 = (al2be2 - gama * gama);
540 sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) * gmx_invsqrt(al2be2*al2be2);
543 /* --- Step4 A3' --- */
544 tmp2 = 1.0 - sinthe * sinthe;
545 costhe = tmp2*gmx_invsqrt(tmp2);
546 xa3d = -ya2d * sinthe;
547 ya3d = ya2d * costhe;
549 xb3d = xb2d * costhe - yb2d * sinthe;
550 yb3d = xb2d * sinthe + yb2d * costhe;
552 xc3d = -xb2d * costhe - yc2d * sinthe;
553 yc3d = -xb2d * sinthe + yc2d * costhe;
557 /* --- Step5 A3 --- */
558 xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
559 ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
560 za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
561 xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
562 yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
563 zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
564 xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
565 yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
566 zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
568 after[ow1] = xcom + xa3;
569 after[ow1 + 1] = ycom + ya3;
570 after[ow1 + 2] = zcom + za3;
571 after[hw2] = xcom + xb3;
572 after[hw2 + 1] = ycom + yb3;
573 after[hw2 + 2] = zcom + zb3;
574 after[hw3] = xcom + xc3;
575 after[hw3 + 1] = ycom + yc3;
576 after[hw3 + 2] = zcom + zc3;
581 rvec_inc(after+hw2, sh_hw2);
582 rvec_inc(after+hw3, sh_hw3);
594 /* 9 flops, counted with the virial */
598 v[ow1] += dax*invdts;
599 v[ow1 + 1] += day*invdts;
600 v[ow1 + 2] += daz*invdts;
601 v[hw2] += dbx*invdts;
602 v[hw2 + 1] += dby*invdts;
603 v[hw2 + 2] += dbz*invdts;
604 v[hw3] += dcx*invdts;
605 v[hw3 + 1] += dcy*invdts;
606 v[hw3 + 2] += dcz*invdts;
610 if (ow1 < CalcVirAtomEnd)
621 vir_r_m_dr[XX][XX] -= b4[ow1 ]*mdax + (b4[ow1 ]+xb0)*mdbx + (b4[ow1 ]+xc0)*mdcx;
622 vir_r_m_dr[XX][YY] -= b4[ow1 ]*mday + (b4[ow1 ]+xb0)*mdby + (b4[ow1 ]+xc0)*mdcy;
623 vir_r_m_dr[XX][ZZ] -= b4[ow1 ]*mdaz + (b4[ow1 ]+xb0)*mdbz + (b4[ow1 ]+xc0)*mdcz;
624 vir_r_m_dr[YY][XX] -= b4[ow1+1]*mdax + (b4[ow1+1]+yb0)*mdbx + (b4[ow1+1]+yc0)*mdcx;
625 vir_r_m_dr[YY][YY] -= b4[ow1+1]*mday + (b4[ow1+1]+yb0)*mdby + (b4[ow1+1]+yc0)*mdcy;
626 vir_r_m_dr[YY][ZZ] -= b4[ow1+1]*mdaz + (b4[ow1+1]+yb0)*mdbz + (b4[ow1+1]+yc0)*mdcz;
627 vir_r_m_dr[ZZ][XX] -= b4[ow1+2]*mdax + (b4[ow1+2]+zb0)*mdbx + (b4[ow1+2]+zc0)*mdcx;
628 vir_r_m_dr[ZZ][YY] -= b4[ow1+2]*mday + (b4[ow1+2]+zb0)*mdby + (b4[ow1+2]+zc0)*mdcy;
629 vir_r_m_dr[ZZ][ZZ] -= b4[ow1+2]*mdaz + (b4[ow1+2]+zb0)*mdbz + (b4[ow1+2]+zc0)*mdcz;
640 check_cons(debug, "settle", after, ow1, hw2, hw3);