src/gromacs/mdlib/csettle.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  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.
  15
  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.
  20  *
  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.
  27  *
  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.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <math.h>
  41 #include <stdio.h>
  42 #include "vec.h"
  43 #include "constr.h"
  44 #include "gmx_fatal.h"
  45 #include "smalloc.h"
  46
  47 typedef struct
  48 {
  49     real   mO;
  50     real   mH;
  51     double wo;
  52     double wh;
  53     double wohh;
  54     real   dOH;
  55     real   dHH;
  56     real   ra;
  57     real   rb;
  58     real   rc;
  59     real   rc2;
  60     /* For projection */
  61     real   imO;
  62     real   imH;
  63     real   invdOH;
  64     real   invdHH;
  65     matrix invmat;
  66 } settleparam_t;
  67
  68 typedef struct gmx_settledata
  69 {
  70     settleparam_t massw;
  71     settleparam_t mass1;
  72 } t_gmx_settledata;
  73
  74
  75 static void init_proj_matrix(settleparam_t *p,
  76                              real invmO,real invmH,real dOH,real dHH)
  77 {
  78     real   imOn,imHn;
  79     matrix mat;
  80
  81     p->imO = invmO;
  82     p->imH = invmH;
  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.
  86      */
  87     imOn = 1;
  88     imHn = p->imH/p->imO;
  89
  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];
 100
 101     m_inv(mat,p->invmat);
 102
 103     msmul(p->invmat,1/p->imO,p->invmat);
 104
 105     p->invdOH = 1/dOH;
 106     p->invdHH = 1/dHH;
 107 }
 108
 109 static void settleparam_init(settleparam_t *p,
 110                              real mO,real mH,real invmO,real invmH,
 111                              real dOH,real dHH)
 112 {
 113     p->mO     = mO;
 114     p->mH     = mH;
 115     p->wo     = mO;
 116     p->wh     = mH;
 117     p->wohh   = mO + 2.0*mH;
 118     p->dOH    = dOH;
 119     p->dHH    = dHH;
 120     p->rc     = dHH/2.0;
 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;
 123     p->rc2    = dHH;
 124
 125     p->wo    /= p->wohh;
 126     p->wh    /= p->wohh;
 127
 128     /* For projection: connection matrix inversion */
 129     init_proj_matrix(p,invmO,invmH,dOH,dHH);
 130
 131     if (debug)
 132     {
 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);
 137     }
 138 }
 139
 140 gmx_settledata_t settle_init(real mO,real mH,real invmO,real invmH,
 141                              real dOH,real dHH)
 142 {
 143     gmx_settledata_t settled;
 144
 145     snew(settled,1);
 146
 147     settleparam_init(&settled->massw,mO,mH,invmO,invmH,dOH,dHH);
 148
 149     settleparam_init(&settled->mass1,1.0,1.0,1.0,1.0,dOH,dHH);
 150
 151     return settled;
 152 }
 153
 154 #ifdef DEBUG
 155 static void check_cons(FILE *fp,char *title,real x[],int OW1,int HW2,int HW3)
 156 {
 157   rvec dOH1,dOH2,dHH;
 158   int  m;
 159
 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];
 164   }
 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));
 167 }
 168 #endif
 169
 170
 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)
 176 {
 177     /* Settle for projection out constraint components
 178      * of derivatives of the coordinates.
 179      * Berk Hess 2008-1-10
 180      */
 181
 182     settleparam_t *p;
 183     real   imO,imH,dOH,dHH,invdOH,invdHH;
 184     matrix invmat;
 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;
 189     real   kfacOH,kfacHH;
 190
 191     if (econq == econqForce)
 192     {
 193         p = &settled->mass1;
 194     }
 195     else
 196     {
 197         p = &settled->massw;
 198     }
 199     imO    = p->imO;
 200     imH    = p->imH;
 201     copy_mat(p->invmat,invmat);
 202     dOH    = p->dOH;
 203     dHH    = p->dHH;
 204     invdOH = p->invdOH;
 205     invdHH = p->invdHH;
 206
 207     veta = vetavar->veta;
 208     vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
 209
 210 #ifdef PRAGMAS
 211 #pragma ivdep
 212 #endif
 213
 214     for (i=0; i<nsettle; i++)
 215     {
 216         ow1 = iatoms[i*4+1];
 217         hw2 = iatoms[i*4+2];
 218         hw3 = iatoms[i*4+3];
 219
 220
 221         for(m=0; m<DIM; m++)
 222         {
 223             /* in the velocity case, these are the velocities, so we
 224                need to modify with the pressure control velocities! */
 225
 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];
 229
 230         }
 231         /* 27 flops */
 232
 233         for(m=0; m<DIM; m++)
 234         {
 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;
 238         }
 239         /* 18 flops */
 240
 241         /* Determine the projections of der(modified) on the bonds */
 242         clear_rvec(dc);
 243         for(m=0; m<DIM; m++)
 244         {
 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];
 248         }
 249         /* 27 flops */
 250
 251         /* Determine the correction for the three bonds */
 252         mvmul(invmat,dc,fc);
 253         /* 15 flops */
 254
 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);
 258         /* 7? flops */
 259
 260         /* Subtract the corrections from derp */
 261         for(m=0; m<DIM; m++)
 262         {
 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]);
 266         }
 267
 268         /* 45 flops */
 269
 270         if (bCalcVir)
 271         {
 272             /* Determining r \dot m der is easy,
 273              * since fc contains the mass weighted corrections for der.
 274              */
 275
 276             for(m=0; m<DIM; m++)
 277             {
 278                 for(m2=0; m2<DIM; m2++)
 279                 {
 280                     rmdder[m][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];
 284                 }
 285             }
 286         }
 287     }
 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);
 291 }
 292
 293
 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)
 296 {
 297   real bondsq[3];
 298   real bond[9];
 299   real invmass[3];
 300   real M2[3];
 301   int iconv;
 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;
 306
 307   invmass[0]=1.0/mO;
 308   invmass[1]=1.0/mH;
 309   invmass[2]=1.0/mH;
 310
 311   bondsq[0]=dOH*dOH;
 312   bondsq[1]=bondsq[0];
 313   bondsq[2]=dHH*dHH;
 314
 315   M2[0]=1.0/(2.0*(invmass[0]+invmass[1]));
 316   M2[1]=M2[0];
 317   M2[2]=1.0/(2.0*(invmass[1]+invmass[2]));
 318
 319   for(ll=0;ll<3;ll++) {
 320     l3=3*ll;
 321     ix=3*iatom[ll];
 322     jx=3*jatom[ll];
 323     for(i=0;i<3;i++)
 324       bond[l3+i]= b4[ix+i] - b4[jx+i];
 325   }
 326
 327   for(i=0,iconv=0;i<1000 && iconv<3; i++) {
 328     for(ll=0;ll<3;ll++) {
 329       ii = iatom[ll];
 330       jj = jatom[ll];
 331       l3 = 3*ll;
 332       ix = 3*ii;
 333       jx = 3*jj;
 334       iy = ix+1;
 335       jy = jx+1;
 336       iz = ix+2;
 337       jz = jx+2;
 338
 339       rijx = bond[l3];
 340       rijy = bond[l3+1];
 341       rijz = bond[l3+2];
 342
 343
 344       tx   = after[ix]-after[jx];
 345       ty   = after[iy]-after[jy];
 346       tz   = after[iz]-after[jz];
 347
 348       rp   = tx*tx+ty*ty+tz*tz;
 349       diff = bondsq[ll] - rp;
 350
 351       if(fabs(diff)<1e-8) {
 352         iconv++;
 353       } else {
 354         rp = rijx*tx+rijy*ty+rijz*tz;
 355         if(rp<1e-8) {
 356           return -1;
 357         }
 358         acor = diff*M2[ll]/rp;
 359         im           = invmass[ii];
 360         jm           = invmass[jj];
 361         tx           = rijx*acor;
 362         ty           = rijy*acor;
 363         tz           = rijz*acor;
 364         after[ix] += tx*im;
 365         after[iy] += ty*im;
 366         after[iz] += tz*im;
 367         after[jx] -= tx*jm;
 368         after[jy] -= ty*jm;
 369         after[jz] -= tz*jm;
 370       }
 371     }
 372   }
 373   return 0;
 374 }
 375
 376
 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)
 380 {
 381     /* ***************************************************************** */
 382     /*                                                               ** */
 383     /*    Subroutine : setlep - reset positions of TIP3P waters      ** */
 384     /*    Author : Shuichi Miyamoto                                  ** */
 385     /*    Date of last update : Oct. 1, 1992                         ** */
 386     /*                                                               ** */
 387     /*    Reference for the SETTLE algorithm                         ** */
 388     /*           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
 389     /*                                                               ** */
 390     /* ***************************************************************** */
 391
 392     /* Initialized data */
 393     settleparam_t *p;
 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. */
 397     double wo,wh,wohh;
 398     real   ra,rb,rc,rc2,dOH,dHH;
 399
 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;
 411     real t1,t2;
 412     real dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
 413     real mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
 414
 415     int doshake;
 416
 417     int i, shakeret, ow1, hw2, hw3;
 418
 419     *error = -1;
 420
 421     p = &settled->massw;
 422     mO   = p->mO;
 423     mH   = p->mH;
 424     wo   = p->wo;
 425     wh   = p->wh;
 426     wohh = p->wohh;
 427     rc   = p->rc;
 428     ra   = p->ra;
 429     rb   = p->rb;
 430     rc2  = p->rc2;
 431     dOH  = p->dOH;
 432     dHH  = p->dHH;
 433
 434     mOs  = mO / vetavar->rvscale;
 435     mHs  = mH / vetavar->rvscale;
 436     invdts = invdt/(vetavar->rscale);
 437
 438 #ifdef PRAGMAS
 439 #pragma ivdep
 440 #endif
 441   for (i = 0; i < nsettle; ++i) {
 442     doshake = 0;
 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];
 453     /* 6 flops */
 454
 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);
 458     /* 12 flops */
 459
 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;
 469     /* 9 flops */
 470
 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;
 480     /* 27 flops */
 481
 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);
 485
 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;
 495     /* 24 flops */
 496
 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;
 501     /*
 502     xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
 503     ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
 504     */
 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;
 512     /* 65 flops */
 513
 514     sinphi = za1d / ra;
 515     tmp    = 1.0 - sinphi * sinphi;
 516     if (tmp <= 0) {
 517       *error = i;
 518       doshake = 1;
 519       cosphi = 0;
 520     }
 521     else
 522       cosphi = tmp*gmx_invsqrt(tmp);
 523     sinpsi = (zb1d - zc1d) / (rc2 * cosphi);
 524     tmp2   = 1.0 - sinpsi * sinpsi;
 525     if (tmp2 <= 0) {
 526       *error = i;
 527       doshake = 1;
 528       cospsi = 0;
 529     }
 530     else
 531       cospsi = tmp2*gmx_invsqrt(tmp2);
 532     /* 46 flops */
 533
 534     if(!doshake) {
 535       ya2d =  ra * cosphi;
 536       xb2d = -rc * cospsi;
 537       t1   = -rb * cosphi;
 538       t2   =  rc * sinpsi * sinphi;
 539       yb2d =  t1 - t2;
 540       yc2d =  t1 + t2;
 541       /* 7 flops */
 542
 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;
 550       /* 47 flops */
 551
 552       /*  --- Step4  A3' --- */
 553       tmp2  = 1.0 - sinthe *sinthe;
 554       costhe = tmp2*gmx_invsqrt(tmp2);
 555       xa3d = -ya2d * sinthe;
 556       ya3d = ya2d * costhe;
 557       za3d = za1d;
 558       xb3d = xb2d * costhe - yb2d * sinthe;
 559       yb3d = xb2d * sinthe + yb2d * costhe;
 560       zb3d = zb1d;
 561       xc3d = -xb2d * costhe - yc2d * sinthe;
 562       yc3d = -xb2d * sinthe + yc2d * costhe;
 563       zc3d = zc1d;
 564       /* 26 flops */
 565
 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;
 576       /* 45 flops */
 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;
 586       /* 9 flops */
 587
 588       dax = xa3 - xa1;
 589       day = ya3 - ya1;
 590       daz = za3 - za1;
 591       dbx = xb3 - xb1;
 592       dby = yb3 - yb1;
 593       dbz = zb3 - zb1;
 594       dcx = xc3 - xc1;
 595       dcy = yc3 - yc1;
 596       dcz = zc3 - zc1;
 597       /* 9 flops, counted with the virial */
 598
 599       if (v) {
 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;
 609           /* 3*6 flops */
 610       }
 611
 612       if (bCalcVir) {
 613           mdax = mOs*dax;
 614           mday = mOs*day;
 615           mdaz = mOs*daz;
 616           mdbx = mHs*dbx;
 617           mdby = mHs*dby;
 618           mdbz = mHs*dbz;
 619           mdcx = mHs*dcx;
 620           mdcy = mHs*dcy;
 621           mdcz = mHs*dcz;
 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;
 631         /* 3*24 - 9 flops */
 632       }
 633     } else {
 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)
 637         *error=i;
 638     }
 639 #ifdef DEBUG
 640     if (debug)
 641     {
 642         check_cons(debug,"settle",after,ow1,hw2,hw3);
 643     }
 644 #endif
 645   }
 646 }