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"
68 typedef struct gmx_settledata
75 static void init_proj_matrix(settleparam_t *p,
76 real invmO,real invmH,real dOH,real dHH)
83 /* We normalize the inverse masses with imO for the matrix inversion.
84 * so we can keep using masses of almost zero for frozen particles,
85 * without running out of the float range in m_inv.
90 /* Construct the constraint coupling matrix */
91 mat[0][0] = imOn + imHn;
92 mat[0][1] = imOn*(1 - 0.5*dHH*dHH/(dOH*dOH));
93 mat[0][2] = imHn*0.5*dHH/dOH;
94 mat[1][1] = mat[0][0];
95 mat[1][2] = mat[0][2];
96 mat[2][2] = imHn + imHn;
97 mat[1][0] = mat[0][1];
98 mat[2][0] = mat[0][2];
99 mat[2][1] = mat[1][2];
101 m_inv(mat,p->invmat);
103 msmul(p->invmat,1/p->imO,p->invmat);
109 static void settleparam_init(settleparam_t *p,
110 real mO,real mH,real invmO,real invmH,
117 p->wohh = mO + 2.0*mH;
121 p->ra = 2.0*p->wh*sqrt(dOH*dOH - p->rc*p->rc)/p->wohh;
122 p->rb = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
128 /* For projection: connection matrix inversion */
129 init_proj_matrix(p,invmO,invmH,dOH,dHH);
133 fprintf(debug,"wo = %g, wh =%g, wohh = %g, rc = %g, ra = %g\n",
134 p->wo,p->wh,p->wohh,p->rc,p->ra);
135 fprintf(debug,"rb = %g, rc2 = %g, dHH = %g, dOH = %g\n",
136 p->rb,p->rc2,p->dHH,p->dOH);
140 gmx_settledata_t settle_init(real mO,real mH,real invmO,real invmH,
143 gmx_settledata_t settled;
147 settleparam_init(&settled->massw,mO,mH,invmO,invmH,dOH,dHH);
149 settleparam_init(&settled->mass1,1.0,1.0,1.0,1.0,dOH,dHH);
155 static void check_cons(FILE *fp,char *title,real x[],int OW1,int HW2,int HW3)
160 for(m=0; (m<DIM); m++) {
161 dOH1[m]=x[OW1+m]-x[HW2+m];
162 dOH2[m]=x[OW1+m]-x[HW3+m];
163 dHH[m]=x[HW2+m]-x[HW3+m];
165 fprintf(fp,"%10s, OW1=%3d, HW2=%3d, HW3=%3d, dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
166 title,OW1/DIM,HW2/DIM,HW3/DIM,norm(dOH1),norm(dOH2),norm(dHH));
171 void settle_proj(FILE *fp,
172 gmx_settledata_t settled,int econq,
173 int nsettle, t_iatom iatoms[],rvec x[],
174 rvec *der,rvec *derp,
175 gmx_bool bCalcVir,tensor rmdder,t_vetavars *vetavar)
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 if (econq == econqForce)
201 copy_mat(p->invmat,invmat);
207 veta = vetavar->veta;
208 vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
214 for (i=0; i<nsettle; i++)
223 /* in the velocity case, these are the velocities, so we
224 need to modify with the pressure control velocities! */
226 derm[0][m] = vscale_nhc*der[ow1][m] + veta*x[ow1][m];
227 derm[1][m] = vscale_nhc*der[hw2][m] + veta*x[hw2][m];
228 derm[2][m] = vscale_nhc*der[hw3][m] + veta*x[hw3][m];
235 roh2[m] = (x[ow1][m] - x[hw2][m])*invdOH;
236 roh3[m] = (x[ow1][m] - x[hw3][m])*invdOH;
237 rhh [m] = (x[hw2][m] - x[hw3][m])*invdHH;
241 /* Determine the projections of der(modified) on the bonds */
245 dc[0] += (derm[0][m] - derm[1][m])*roh2[m];
246 dc[1] += (derm[0][m] - derm[2][m])*roh3[m];
247 dc[2] += (derm[1][m] - derm[2][m])*rhh [m];
251 /* Determine the correction for the three bonds */
255 /* divide velocity by vscale_nhc for determining constrained velocities, since they
256 have not yet been multiplied */
257 svmul(1.0/vscale_nhc,fc,fcv);
260 /* Subtract the corrections from derp */
263 derp[ow1][m] -= imO*( fcv[0]*roh2[m] + fcv[1]*roh3[m]);
264 derp[hw2][m] -= imH*(-fcv[0]*roh2[m] + fcv[2]*rhh [m]);
265 derp[hw3][m] -= imH*(-fcv[1]*roh3[m] - fcv[2]*rhh [m]);
272 /* Determining r \dot m der is easy,
273 * since fc contains the mass weighted corrections for der.
278 for(m2=0; m2<DIM; m2++)
281 dOH*roh2[m]*roh2[m2]*fcv[0] +
282 dOH*roh3[m]*roh3[m2]*fcv[1] +
283 dHH*rhh [m]*rhh [m2]*fcv[2];
288 /* conrect rmdder, which will be used to calcualate the virial; we need to use
289 the unscaled multipliers in the virial */
290 msmul(rmdder,1.0/vetavar->vscale,rmdder);
294 /* Our local shake routine to be used when settle breaks down due to a zero determinant */
295 static int xshake(real b4[], real after[], real dOH, real dHH, real mO, real mH)
302 int iatom[3]={0,0,1};
303 int jatom[3]={1,2,2};
304 real rijx,rijy,rijz,tx,ty,tz,im,jm,acor,rp,diff;
305 int i,ll,ii,jj,l3,ix,iy,iz,jx,jy,jz,conv;
315 M2[0]=1.0/(2.0*(invmass[0]+invmass[1]));
317 M2[2]=1.0/(2.0*(invmass[1]+invmass[2]));
319 for(ll=0;ll<3;ll++) {
324 bond[l3+i]= b4[ix+i] - b4[jx+i];
327 for(i=0,iconv=0;i<1000 && iconv<3; i++) {
328 for(ll=0;ll<3;ll++) {
344 tx = after[ix]-after[jx];
345 ty = after[iy]-after[jy];
346 tz = after[iz]-after[jz];
348 rp = tx*tx+ty*ty+tz*tz;
349 diff = bondsq[ll] - rp;
351 if(fabs(diff)<1e-8) {
354 rp = rijx*tx+rijy*ty+rijz*tz;
358 acor = diff*M2[ll]/rp;
377 void csettle(gmx_settledata_t settled,
378 int nsettle, t_iatom iatoms[],real b4[], real after[],
379 real invdt,real *v,gmx_bool bCalcVir,tensor rmdr,int *error,t_vetavars *vetavar)
381 /* ***************************************************************** */
383 /* Subroutine : setlep - reset positions of TIP3P waters ** */
384 /* Author : Shuichi Miyamoto ** */
385 /* Date of last update : Oct. 1, 1992 ** */
387 /* Reference for the SETTLE algorithm ** */
388 /* S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
390 /* ***************************************************************** */
392 /* Initialized data */
394 real mO,mH,mOs,mHs,invdts;
395 /* These three weights need have double precision. Using single precision
396 * can result in huge velocity and pressure deviations. */
398 real ra,rb,rc,rc2,dOH,dHH;
400 /* Local variables */
401 real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2;
402 real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22,
403 trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe,
404 cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd,
405 xb0, yb0, zb0, xc0, yc0, zc0, xa1;
406 real ya1, za1, xb1, yb1;
407 real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3,
408 xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d,
409 za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d,
410 xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d;
412 real dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
413 real mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
417 int i, shakeret, ow1, hw2, hw3;
434 mOs = mO / vetavar->rvscale;
435 mHs = mH / vetavar->rvscale;
436 invdts = invdt/(vetavar->rscale);
441 for (i = 0; i < nsettle; ++i) {
443 /* --- Step1 A1' --- */
444 ow1 = iatoms[i*4+1] * 3;
445 hw2 = iatoms[i*4+2] * 3;
446 hw3 = iatoms[i*4+3] * 3;
447 xb0 = b4[hw2 ] - b4[ow1];
448 yb0 = b4[hw2 + 1] - b4[ow1 + 1];
449 zb0 = b4[hw2 + 2] - b4[ow1 + 2];
450 xc0 = b4[hw3 ] - b4[ow1];
451 yc0 = b4[hw3 + 1] - b4[ow1 + 1];
452 zc0 = b4[hw3 + 2] - b4[ow1 + 2];
455 xcom = (after[ow1 ] * wo + (after[hw2 ] + after[hw3 ]) * wh);
456 ycom = (after[ow1 + 1] * wo + (after[hw2 + 1] + after[hw3 + 1]) * wh);
457 zcom = (after[ow1 + 2] * wo + (after[hw2 + 2] + after[hw3 + 2]) * wh);
460 xa1 = after[ow1 ] - xcom;
461 ya1 = after[ow1 + 1] - ycom;
462 za1 = after[ow1 + 2] - zcom;
463 xb1 = after[hw2 ] - xcom;
464 yb1 = after[hw2 + 1] - ycom;
465 zb1 = after[hw2 + 2] - zcom;
466 xc1 = after[hw3 ] - xcom;
467 yc1 = after[hw3 + 1] - ycom;
468 zc1 = after[hw3 + 2] - zcom;
471 xakszd = yb0 * zc0 - zb0 * yc0;
472 yakszd = zb0 * xc0 - xb0 * zc0;
473 zakszd = xb0 * yc0 - yb0 * xc0;
474 xaksxd = ya1 * zakszd - za1 * yakszd;
475 yaksxd = za1 * xakszd - xa1 * zakszd;
476 zaksxd = xa1 * yakszd - ya1 * xakszd;
477 xaksyd = yakszd * zaksxd - zakszd * yaksxd;
478 yaksyd = zakszd * xaksxd - xakszd * zaksxd;
479 zaksyd = xakszd * yaksxd - yakszd * xaksxd;
482 axlng = gmx_invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
483 aylng = gmx_invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
484 azlng = gmx_invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
486 trns11 = xaksxd * axlng;
487 trns21 = yaksxd * axlng;
488 trns31 = zaksxd * axlng;
489 trns12 = xaksyd * aylng;
490 trns22 = yaksyd * aylng;
491 trns32 = zaksyd * aylng;
492 trns13 = xakszd * azlng;
493 trns23 = yakszd * azlng;
494 trns33 = zakszd * azlng;
497 xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
498 yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
499 xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
500 yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
502 xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
503 ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
505 za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
506 xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
507 yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
508 zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
509 xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
510 yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
511 zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
515 tmp = 1.0 - sinphi * sinphi;
522 cosphi = tmp*gmx_invsqrt(tmp);
523 sinpsi = (zb1d - zc1d) / (rc2 * cosphi);
524 tmp2 = 1.0 - sinpsi * sinpsi;
531 cospsi = tmp2*gmx_invsqrt(tmp2);
538 t2 = rc * sinpsi * sinphi;
543 /* --- Step3 al,be,ga --- */
544 alpa = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
545 beta = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
546 gama = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
547 al2be2 = alpa * alpa + beta * beta;
548 tmp2 = (al2be2 - gama * gama);
549 sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) / al2be2;
552 /* --- Step4 A3' --- */
553 tmp2 = 1.0 - sinthe *sinthe;
554 costhe = tmp2*gmx_invsqrt(tmp2);
555 xa3d = -ya2d * sinthe;
556 ya3d = ya2d * costhe;
558 xb3d = xb2d * costhe - yb2d * sinthe;
559 yb3d = xb2d * sinthe + yb2d * costhe;
561 xc3d = -xb2d * costhe - yc2d * sinthe;
562 yc3d = -xb2d * sinthe + yc2d * costhe;
566 /* --- Step5 A3 --- */
567 xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
568 ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
569 za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
570 xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
571 yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
572 zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
573 xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
574 yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
575 zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
577 after[ow1] = xcom + xa3;
578 after[ow1 + 1] = ycom + ya3;
579 after[ow1 + 2] = zcom + za3;
580 after[hw2] = xcom + xb3;
581 after[hw2 + 1] = ycom + yb3;
582 after[hw2 + 2] = zcom + zb3;
583 after[hw3] = xcom + xc3;
584 after[hw3 + 1] = ycom + yc3;
585 after[hw3 + 2] = zcom + zc3;
597 /* 9 flops, counted with the virial */
600 v[ow1] += dax*invdts;
601 v[ow1 + 1] += day*invdts;
602 v[ow1 + 2] += daz*invdts;
603 v[hw2] += dbx*invdts;
604 v[hw2 + 1] += dby*invdts;
605 v[hw2 + 2] += dbz*invdts;
606 v[hw3] += dcx*invdts;
607 v[hw3 + 1] += dcy*invdts;
608 v[hw3 + 2] += dcz*invdts;
622 rmdr[XX][XX] -= b4[ow1]*mdax + b4[hw2]*mdbx + b4[hw3]*mdcx;
623 rmdr[XX][YY] -= b4[ow1]*mday + b4[hw2]*mdby + b4[hw3]*mdcy;
624 rmdr[XX][ZZ] -= b4[ow1]*mdaz + b4[hw2]*mdbz + b4[hw3]*mdcz;
625 rmdr[YY][XX] -= b4[ow1+1]*mdax + b4[hw2+1]*mdbx + b4[hw3+1]*mdcx;
626 rmdr[YY][YY] -= b4[ow1+1]*mday + b4[hw2+1]*mdby + b4[hw3+1]*mdcy;
627 rmdr[YY][ZZ] -= b4[ow1+1]*mdaz + b4[hw2+1]*mdbz + b4[hw3+1]*mdcz;
628 rmdr[ZZ][XX] -= b4[ow1+2]*mdax + b4[hw2+2]*mdbx + b4[hw3+2]*mdcx;
629 rmdr[ZZ][YY] -= b4[ow1+2]*mday + b4[hw2+2]*mdby + b4[hw3+2]*mdcy;
630 rmdr[ZZ][ZZ] -= b4[ow1+2]*mdaz + b4[hw2+2]*mdbz + b4[hw3+2]*mdcz;
634 /* If we couldn't settle this water, try a simplified iterative shake instead */
635 /* no pressure control in here yet */
636 if(xshake(b4+ow1,after+ow1,dOH,dHH,mO,mH)!=0)
642 check_cons(debug,"settle",after,ow1,hw2,hw3);