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",
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)
157 for(m=0; (m<DIM); m++) {
158 dOH1[m]=x[OW1+m]-x[HW2+m];
159 dOH2[m]=x[OW1+m]-x[HW3+m];
160 dHH[m]=x[HW2+m]-x[HW3+m];
162 fprintf(fp,"%10s, OW1=%3d, HW2=%3d, HW3=%3d, dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
163 title,OW1/DIM,HW2/DIM,HW3/DIM,norm(dOH1),norm(dOH2),norm(dHH));
168 void settle_proj(FILE *fp,
169 gmx_settledata_t settled,int econq,
170 int nsettle, t_iatom iatoms[],
173 rvec *der,rvec *derp,
174 int calcvir_atom_end,tensor vir_r_m_dder,
177 /* Settle for projection out constraint components
178 * of derivatives of the coordinates.
179 * Berk Hess 2008-1-10
183 real imO,imH,dOH,dHH,invdOH,invdHH;
185 int i,m,m2,ow1,hw2,hw3;
186 rvec roh2,roh3,rhh,dc,fc,fcv;
187 rvec derm[3],derpm[3];
188 real invvscale,vscale_nhc,veta;
191 calcvir_atom_end *= DIM;
193 if (econq == econqForce)
203 copy_mat(p->invmat,invmat);
209 veta = vetavar->veta;
210 vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
216 for (i=0; i<nsettle; i++)
225 /* in the velocity case, these are the velocities, so we
226 need to modify with the pressure control velocities! */
228 derm[0][m] = vscale_nhc*der[ow1][m] + veta*x[ow1][m];
229 derm[1][m] = vscale_nhc*der[hw2][m] + veta*x[hw2][m];
230 derm[2][m] = vscale_nhc*der[hw3][m] + veta*x[hw3][m];
237 rvec_sub(x[ow1],x[hw2],roh2);
238 rvec_sub(x[ow1],x[hw3],roh3);
239 rvec_sub(x[hw2],x[hw3],rhh);
243 pbc_dx_aiuc(pbc,x[ow1],x[hw2],roh2);
244 pbc_dx_aiuc(pbc,x[ow1],x[hw3],roh3);
245 pbc_dx_aiuc(pbc,x[hw2],x[hw3],rhh);
247 svmul(invdOH,roh2,roh2);
248 svmul(invdOH,roh3,roh3);
249 svmul(invdHH,rhh,rhh);
252 /* Determine the projections of der(modified) on the bonds */
256 dc[0] += (derm[0][m] - derm[1][m])*roh2[m];
257 dc[1] += (derm[0][m] - derm[2][m])*roh3[m];
258 dc[2] += (derm[1][m] - derm[2][m])*rhh [m];
262 /* Determine the correction for the three bonds */
266 /* divide velocity by vscale_nhc for determining constrained velocities, since they
267 have not yet been multiplied */
268 svmul(1.0/vscale_nhc,fc,fcv);
271 /* Subtract the corrections from derp */
274 derp[ow1][m] -= imO*( fcv[0]*roh2[m] + fcv[1]*roh3[m]);
275 derp[hw2][m] -= imH*(-fcv[0]*roh2[m] + fcv[2]*rhh [m]);
276 derp[hw3][m] -= imH*(-fcv[1]*roh3[m] - fcv[2]*rhh [m]);
281 if (ow1 < calcvir_atom_end)
283 /* Determining r \dot m der is easy,
284 * since fc contains the mass weighted corrections for der.
289 for(m2=0; m2<DIM; m2++)
291 vir_r_m_dder[m][m2] +=
292 dOH*roh2[m]*roh2[m2]*fcv[0] +
293 dOH*roh3[m]*roh3[m2]*fcv[1] +
294 dHH*rhh [m]*rhh [m2]*fcv[2];
300 if (calcvir_atom_end > 0)
302 /* Correct r_m_dder, which will be used to calcualate the virial;
303 * we need to use the unscaled multipliers in the virial.
305 msmul(vir_r_m_dder,1.0/vetavar->vscale,vir_r_m_dder);
310 void csettle(gmx_settledata_t settled,
311 int nsettle, t_iatom iatoms[],
313 real b4[], real after[],
314 real invdt,real *v,int CalcVirAtomEnd,
319 /* ***************************************************************** */
321 /* Subroutine : setlep - reset positions of TIP3P waters ** */
322 /* Author : Shuichi Miyamoto ** */
323 /* Date of last update : Oct. 1, 1992 ** */
325 /* Reference for the SETTLE algorithm ** */
326 /* S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
328 /* ***************************************************************** */
330 /* Initialized data */
332 real wh,ra,rb,rc,irc2;
335 /* Local variables */
336 real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2;
337 real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22,
338 trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe,
339 cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd,
340 xb0, yb0, zb0, xc0, yc0, zc0, xa1;
341 real ya1, za1, xb1, yb1;
342 real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3,
343 xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d,
344 za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d,
345 xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d;
347 real dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
348 real mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
352 int i, ow1, hw2, hw3;
354 rvec dx,sh_hw2={0,0,0},sh_hw3={0,0,0};
369 mOs = p->mO / vetavar->rvscale;
370 mHs = p->mH / vetavar->rvscale;
371 invdts = invdt / vetavar->rscale;
376 for (i = 0; i < nsettle; ++i) {
378 /* --- Step1 A1' --- */
379 ow1 = iatoms[i*4+1] * 3;
380 hw2 = iatoms[i*4+2] * 3;
381 hw3 = iatoms[i*4+3] * 3;
384 xb0 = b4[hw2 + XX] - b4[ow1 + XX];
385 yb0 = b4[hw2 + YY] - b4[ow1 + YY];
386 zb0 = b4[hw2 + ZZ] - b4[ow1 + ZZ];
387 xc0 = b4[hw3 + XX] - b4[ow1 + XX];
388 yc0 = b4[hw3 + YY] - b4[ow1 + YY];
389 zc0 = b4[hw3 + ZZ] - b4[ow1 + ZZ];
392 rvec_sub(after+hw2,after+ow1,doh2);
393 rvec_sub(after+hw3,after+ow1,doh3);
398 pbc_dx_aiuc(pbc,b4+hw2,b4+ow1,dx);
402 pbc_dx_aiuc(pbc,b4+hw3,b4+ow1,dx);
407 /* Tedious way of doing pbc */
408 is = pbc_dx_aiuc(pbc,after+hw2,after+ow1,doh2);
415 sh_hw2[XX] = after[hw2 + XX] - (after[ow1 + XX] + doh2[XX]);
416 sh_hw2[YY] = after[hw2 + YY] - (after[ow1 + YY] + doh2[YY]);
417 sh_hw2[ZZ] = after[hw2 + ZZ] - (after[ow1 + ZZ] + doh2[ZZ]);
418 rvec_dec(after+hw2,sh_hw2);
420 is = pbc_dx_aiuc(pbc,after+hw3,after+ow1,doh3);
427 sh_hw3[XX] = after[hw3 + XX] - (after[ow1 + XX] + doh3[XX]);
428 sh_hw3[YY] = after[hw3 + YY] - (after[ow1 + YY] + doh3[YY]);
429 sh_hw3[ZZ] = after[hw3 + ZZ] - (after[ow1 + ZZ] + doh3[ZZ]);
430 rvec_dec(after+hw3,sh_hw3);
434 /* Not calculating the center of mass using the oxygen position
435 * and the O-H distances, as done below, will make SETTLE
436 * the largest source of energy drift for simulations of water,
437 * as then the oxygen coordinate is multiplied by 0.89 at every step,
438 * which can then transfer a systematic rounding to the oxygen velocity.
440 xa1 = -(doh2[XX] + doh3[XX]) * wh;
441 ya1 = -(doh2[YY] + doh3[YY]) * wh;
442 za1 = -(doh2[ZZ] + doh3[ZZ]) * wh;
444 xcom = after[ow1 + XX] - xa1;
445 ycom = after[ow1 + YY] - ya1;
446 zcom = after[ow1 + ZZ] - za1;
448 xb1 = after[hw2 + XX] - xcom;
449 yb1 = after[hw2 + YY] - ycom;
450 zb1 = after[hw2 + ZZ] - zcom;
451 xc1 = after[hw3 + XX] - xcom;
452 yc1 = after[hw3 + YY] - ycom;
453 zc1 = after[hw3 + ZZ] - zcom;
456 xakszd = yb0 * zc0 - zb0 * yc0;
457 yakszd = zb0 * xc0 - xb0 * zc0;
458 zakszd = xb0 * yc0 - yb0 * xc0;
459 xaksxd = ya1 * zakszd - za1 * yakszd;
460 yaksxd = za1 * xakszd - xa1 * zakszd;
461 zaksxd = xa1 * yakszd - ya1 * xakszd;
462 xaksyd = yakszd * zaksxd - zakszd * yaksxd;
463 yaksyd = zakszd * xaksxd - xakszd * zaksxd;
464 zaksyd = xakszd * yaksxd - yakszd * xaksxd;
467 axlng = gmx_invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
468 aylng = gmx_invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
469 azlng = gmx_invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
471 trns11 = xaksxd * axlng;
472 trns21 = yaksxd * axlng;
473 trns31 = zaksxd * axlng;
474 trns12 = xaksyd * aylng;
475 trns22 = yaksyd * aylng;
476 trns32 = zaksyd * aylng;
477 trns13 = xakszd * azlng;
478 trns23 = yakszd * azlng;
479 trns33 = zakszd * azlng;
482 xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
483 yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
484 xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
485 yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
487 xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
488 ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
490 za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
491 xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
492 yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
493 zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
494 xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
495 yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
496 zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
499 sinphi = za1d * gmx_invsqrt(ra*ra);
500 tmp = 1.0 - sinphi * sinphi;
504 tmp2 = gmx_invsqrt(tmp);
506 sinpsi = (zb1d - zc1d) * irc2 * tmp2;
507 tmp2 = 1.0 - sinpsi * sinpsi;
511 cospsi = tmp2*gmx_invsqrt(tmp2);
520 t2 = rc * sinpsi * sinphi;
525 /* --- Step3 al,be,ga --- */
526 alpa = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
527 beta = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
528 gama = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
529 al2be2 = alpa * alpa + beta * beta;
530 tmp2 = (al2be2 - gama * gama);
531 sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) * gmx_invsqrt(al2be2*al2be2);
534 /* --- Step4 A3' --- */
535 tmp2 = 1.0 - sinthe * sinthe;
536 costhe = tmp2*gmx_invsqrt(tmp2);
537 xa3d = -ya2d * sinthe;
538 ya3d = ya2d * costhe;
540 xb3d = xb2d * costhe - yb2d * sinthe;
541 yb3d = xb2d * sinthe + yb2d * costhe;
543 xc3d = -xb2d * costhe - yc2d * sinthe;
544 yc3d = -xb2d * sinthe + yc2d * costhe;
548 /* --- Step5 A3 --- */
549 xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
550 ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
551 za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
552 xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
553 yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
554 zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
555 xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
556 yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
557 zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
559 after[ow1] = xcom + xa3;
560 after[ow1 + 1] = ycom + ya3;
561 after[ow1 + 2] = zcom + za3;
562 after[hw2] = xcom + xb3;
563 after[hw2 + 1] = ycom + yb3;
564 after[hw2 + 2] = zcom + zb3;
565 after[hw3] = xcom + xc3;
566 after[hw3 + 1] = ycom + yc3;
567 after[hw3 + 2] = zcom + zc3;
572 rvec_inc(after+hw2,sh_hw2);
573 rvec_inc(after+hw3,sh_hw3);
585 /* 9 flops, counted with the virial */
588 v[ow1] += dax*invdts;
589 v[ow1 + 1] += day*invdts;
590 v[ow1 + 2] += daz*invdts;
591 v[hw2] += dbx*invdts;
592 v[hw2 + 1] += dby*invdts;
593 v[hw2 + 2] += dbz*invdts;
594 v[hw3] += dcx*invdts;
595 v[hw3 + 1] += dcy*invdts;
596 v[hw3 + 2] += dcz*invdts;
600 if (ow1 < CalcVirAtomEnd) {
610 vir_r_m_dr[XX][XX] -= b4[ow1 ]*mdax + (b4[ow1 ]+xb0)*mdbx + (b4[ow1 ]+xc0)*mdcx;
611 vir_r_m_dr[XX][YY] -= b4[ow1 ]*mday + (b4[ow1 ]+xb0)*mdby + (b4[ow1 ]+xc0)*mdcy;
612 vir_r_m_dr[XX][ZZ] -= b4[ow1 ]*mdaz + (b4[ow1 ]+xb0)*mdbz + (b4[ow1 ]+xc0)*mdcz;
613 vir_r_m_dr[YY][XX] -= b4[ow1+1]*mdax + (b4[ow1+1]+yb0)*mdbx + (b4[ow1+1]+yc0)*mdcx;
614 vir_r_m_dr[YY][YY] -= b4[ow1+1]*mday + (b4[ow1+1]+yb0)*mdby + (b4[ow1+1]+yc0)*mdcy;
615 vir_r_m_dr[YY][ZZ] -= b4[ow1+1]*mdaz + (b4[ow1+1]+yb0)*mdbz + (b4[ow1+1]+yc0)*mdcz;
616 vir_r_m_dr[ZZ][XX] -= b4[ow1+2]*mdax + (b4[ow1+2]+zb0)*mdbx + (b4[ow1+2]+zc0)*mdcx;
617 vir_r_m_dr[ZZ][YY] -= b4[ow1+2]*mday + (b4[ow1+2]+zb0)*mdby + (b4[ow1+2]+zc0)*mdcy;
618 vir_r_m_dr[ZZ][ZZ] -= b4[ow1+2]*mdaz + (b4[ow1+2]+zb0)*mdbz + (b4[ow1+2]+zc0)*mdcz;
627 check_cons(debug,"settle",after,ow1,hw2,hw3);