src/mdlib/csettle.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include <stdio.h>
  44 #include "vec.h"
  45 #include "constr.h"
  46 #include "gmx_fatal.h"
  47 #include "smalloc.h"
  48 #include "pbc.h"
  49
  50 typedef struct
  51 {
  52     real   mO;
  53     real   mH;
  54     real   wh;
  55     real   dOH;
  56     real   dHH;
  57     real   ra;
  58     real   rb;
  59     real   rc;
  60     real   irc2;
  61     /* For projection */
  62     real   imO;
  63     real   imH;
  64     real   invdOH;
  65     real   invdHH;
  66     matrix invmat;
  67 } settleparam_t;
  68
  69 typedef struct gmx_settledata
  70 {
  71     settleparam_t massw;
  72     settleparam_t mass1;
  73 } t_gmx_settledata;
  74
  75
  76 static void init_proj_matrix(settleparam_t *p,
  77                              real invmO,real invmH,real dOH,real dHH)
  78 {
  79     real   imOn,imHn;
  80     matrix mat;
  81
  82     p->imO = invmO;
  83     p->imH = invmH;
  84     /* We normalize the inverse masses with imO for the matrix inversion.
  85      * so we can keep using masses of almost zero for frozen particles,
  86      * without running out of the float range in m_inv.
  87      */
  88     imOn = 1;
  89     imHn = p->imH/p->imO;
  90
  91     /* Construct the constraint coupling matrix */
  92     mat[0][0] = imOn + imHn;
  93     mat[0][1] = imOn*(1 - 0.5*dHH*dHH/(dOH*dOH));
  94     mat[0][2] = imHn*0.5*dHH/dOH;
  95     mat[1][1] = mat[0][0];
  96     mat[1][2] = mat[0][2];
  97     mat[2][2] = imHn + imHn;
  98     mat[1][0] = mat[0][1];
  99     mat[2][0] = mat[0][2];
 100     mat[2][1] = mat[1][2];
 101
 102     m_inv(mat,p->invmat);
 103
 104     msmul(p->invmat,1/p->imO,p->invmat);
 105
 106     p->invdOH = 1/dOH;
 107     p->invdHH = 1/dHH;
 108 }
 109
 110 static void settleparam_init(settleparam_t *p,
 111                              real mO,real mH,real invmO,real invmH,
 112                              real dOH,real dHH)
 113 {
 114     double wohh;
 115
 116     p->mO     = mO;
 117     p->mH     = mH;
 118     wohh      = mO + 2.0*mH;
 119     p->wh     = mH/wohh;
 120     p->dOH    = dOH;
 121     p->dHH    = dHH;
 122     p->rc     = dHH/2.0;
 123     p->ra     = 2.0*mH*sqrt(dOH*dOH - p->rc*p->rc)/wohh;
 124     p->rb     = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
 125     p->irc2   = 1.0/dHH;
 126
 127     /* For projection: connection matrix inversion */
 128     init_proj_matrix(p,invmO,invmH,dOH,dHH);
 129
 130     if (debug)
 131     {
 132         fprintf(debug,"wh =%g, rc = %g, ra = %g\n",
 133                 p->wh,p->rc,p->ra);
 134         fprintf(debug,"rb = %g, irc2 = %g, dHH = %g, dOH = %g\n",
 135                 p->rb,p->irc2,p->dHH,p->dOH);
 136     }
 137 }
 138
 139 gmx_settledata_t settle_init(real mO,real mH,real invmO,real invmH,
 140                              real dOH,real dHH)
 141 {
 142     gmx_settledata_t settled;
 143
 144     snew(settled,1);
 145
 146     settleparam_init(&settled->massw,mO,mH,invmO,invmH,dOH,dHH);
 147
 148     settleparam_init(&settled->mass1,1.0,1.0,1.0,1.0,dOH,dHH);
 149
 150     return settled;
 151 }
 152
 153 #ifdef DEBUG
 154 static void check_cons(FILE *fp,char *title,real x[],int OW1,int HW2,int HW3)
 155 {
 156   rvec dOH1,dOH2,dHH;
 157   int  m;
 158
 159   for(m=0; (m<DIM); m++) {
 160     dOH1[m]=x[OW1+m]-x[HW2+m];
 161     dOH2[m]=x[OW1+m]-x[HW3+m];
 162     dHH[m]=x[HW2+m]-x[HW3+m];
 163   }
 164   fprintf(fp,"%10s, OW1=%3d, HW2=%3d, HW3=%3d,  dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
 165           title,OW1/DIM,HW2/DIM,HW3/DIM,norm(dOH1),norm(dOH2),norm(dHH));
 166 }
 167 #endif
 168
 169
 170 void settle_proj(FILE *fp,
 171                  gmx_settledata_t settled,int econq,
 172                  int nsettle, t_iatom iatoms[],
 173                  const t_pbc *pbc,
 174                  rvec x[],
 175                  rvec *der,rvec *derp,
 176                  int calcvir_atom_end,tensor vir_r_m_dder,
 177                  t_vetavars *vetavar)
 178 {
 179     /* Settle for projection out constraint components
 180      * of derivatives of the coordinates.
 181      * Berk Hess 2008-1-10
 182      */
 183
 184     settleparam_t *p;
 185     real   imO,imH,dOH,dHH,invdOH,invdHH;
 186     matrix invmat;
 187     int    i,m,m2,ow1,hw2,hw3;
 188     rvec   roh2,roh3,rhh,dc,fc,fcv;
 189     rvec   derm[3],derpm[3];
 190     real   invvscale,vscale_nhc,veta;
 191     real   kfacOH,kfacHH;
 192
 193     calcvir_atom_end *= DIM;
 194
 195     if (econq == econqForce)
 196     {
 197         p = &settled->mass1;
 198     }
 199     else
 200     {
 201         p = &settled->massw;
 202     }
 203     imO    = p->imO;
 204     imH    = p->imH;
 205     copy_mat(p->invmat,invmat);
 206     dOH    = p->dOH;
 207     dHH    = p->dHH;
 208     invdOH = p->invdOH;
 209     invdHH = p->invdHH;
 210
 211     veta = vetavar->veta;
 212     vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
 213
 214 #ifdef PRAGMAS
 215 #pragma ivdep
 216 #endif
 217
 218     for (i=0; i<nsettle; i++)
 219     {
 220         ow1 = iatoms[i*4+1];
 221         hw2 = iatoms[i*4+2];
 222         hw3 = iatoms[i*4+3];
 223
 224
 225         for(m=0; m<DIM; m++)
 226         {
 227             /* in the velocity case, these are the velocities, so we
 228                need to modify with the pressure control velocities! */
 229
 230             derm[0][m]  = vscale_nhc*der[ow1][m] + veta*x[ow1][m];
 231             derm[1][m]  = vscale_nhc*der[hw2][m] + veta*x[hw2][m];
 232             derm[2][m]  = vscale_nhc*der[hw3][m] + veta*x[hw3][m];
 233
 234         }
 235         /* 27 flops */
 236
 237         if (pbc == NULL)
 238         {
 239             rvec_sub(x[ow1],x[hw2],roh2);
 240             rvec_sub(x[ow1],x[hw3],roh3);
 241             rvec_sub(x[hw2],x[hw3],rhh);
 242         }
 243         else
 244         {
 245             pbc_dx_aiuc(pbc,x[ow1],x[hw2],roh2);
 246             pbc_dx_aiuc(pbc,x[ow1],x[hw3],roh3);
 247             pbc_dx_aiuc(pbc,x[hw2],x[hw3],rhh);
 248         }
 249         svmul(invdOH,roh2,roh2);
 250         svmul(invdOH,roh3,roh3);
 251         svmul(invdHH,rhh,rhh);
 252         /* 18 flops */
 253
 254         /* Determine the projections of der(modified) on the bonds */
 255         clear_rvec(dc);
 256         for(m=0; m<DIM; m++)
 257         {
 258             dc[0] += (derm[0][m] - derm[1][m])*roh2[m];
 259             dc[1] += (derm[0][m] - derm[2][m])*roh3[m];
 260             dc[2] += (derm[1][m] - derm[2][m])*rhh [m];
 261         }
 262         /* 27 flops */
 263
 264         /* Determine the correction for the three bonds */
 265         mvmul(invmat,dc,fc);
 266         /* 15 flops */
 267
 268         /* divide velocity by vscale_nhc for determining constrained velocities, since they
 269            have not yet been multiplied */
 270         svmul(1.0/vscale_nhc,fc,fcv);
 271         /* 7? flops */
 272
 273         /* Subtract the corrections from derp */
 274         for(m=0; m<DIM; m++)
 275         {
 276             derp[ow1][m] -= imO*( fcv[0]*roh2[m] + fcv[1]*roh3[m]);
 277             derp[hw2][m] -= imH*(-fcv[0]*roh2[m] + fcv[2]*rhh [m]);
 278             derp[hw3][m] -= imH*(-fcv[1]*roh3[m] - fcv[2]*rhh [m]);
 279         }
 280
 281         /* 45 flops */
 282
 283         if (ow1 < calcvir_atom_end)
 284         {
 285             /* Determining r \dot m der is easy,
 286              * since fc contains the mass weighted corrections for der.
 287              */
 288
 289             for(m=0; m<DIM; m++)
 290             {
 291                 for(m2=0; m2<DIM; m2++)
 292                 {
 293                     vir_r_m_dder[m][m2] +=
 294                         dOH*roh2[m]*roh2[m2]*fcv[0] +
 295                         dOH*roh3[m]*roh3[m2]*fcv[1] +
 296                         dHH*rhh [m]*rhh [m2]*fcv[2];
 297                 }
 298             }
 299         }
 300     }
 301
 302     if (calcvir_atom_end > 0)
 303     {
 304         /* Correct r_m_dder, which will be used to calcualate the virial;
 305          * we need to use the unscaled multipliers in the virial.
 306          */
 307         msmul(vir_r_m_dder,1.0/vetavar->vscale,vir_r_m_dder);
 308     }
 309 }
 310
 311
 312 void csettle(gmx_settledata_t settled,
 313              int nsettle, t_iatom iatoms[],
 314              const t_pbc *pbc,
 315              real b4[], real after[],
 316              real invdt,real *v,int CalcVirAtomEnd,
 317              tensor vir_r_m_dr,
 318              int *error,
 319              t_vetavars *vetavar)
 320 {
 321     /* ***************************************************************** */
 322     /*                                                               ** */
 323     /*    Subroutine : setlep - reset positions of TIP3P waters      ** */
 324     /*    Author : Shuichi Miyamoto                                  ** */
 325     /*    Date of last update : Oct. 1, 1992                         ** */
 326     /*                                                               ** */
 327     /*    Reference for the SETTLE algorithm                         ** */
 328     /*           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
 329     /*                                                               ** */
 330     /* ***************************************************************** */
 331
 332     /* Initialized data */
 333     settleparam_t *p;
 334     real   wh,ra,rb,rc,irc2;
 335     real   mOs,mHs,invdts;
 336
 337     /* Local variables */
 338     real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2;
 339     real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22,
 340         trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe,
 341         cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd,
 342         xb0, yb0, zb0, xc0, yc0, zc0, xa1;
 343     real ya1, za1, xb1, yb1;
 344     real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3,
 345         xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d,
 346         za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d,
 347         xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d;
 348     real t1,t2;
 349     real dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
 350     real mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
 351
 352     gmx_bool bOK;
 353
 354     int i, ow1, hw2, hw3;
 355
 356     rvec dx,sh_hw2={0,0,0},sh_hw3={0,0,0};
 357     rvec doh2,doh3;
 358     int  is;
 359
 360     *error = -1;
 361
 362     CalcVirAtomEnd *= 3;
 363
 364     p = &settled->massw;
 365     wh   = p->wh;
 366     rc   = p->rc;
 367     ra   = p->ra;
 368     rb   = p->rb;
 369     irc2 = p->irc2;
 370
 371     mOs    = p->mO / vetavar->rvscale;
 372     mHs    = p->mH / vetavar->rvscale;
 373     invdts = invdt / vetavar->rscale;
 374
 375 #ifdef PRAGMAS
 376 #pragma ivdep
 377 #endif
 378   for (i = 0; i < nsettle; ++i) {
 379     bOK = TRUE;
 380     /*    --- Step1  A1' ---      */
 381     ow1 = iatoms[i*4+1] * 3;
 382     hw2 = iatoms[i*4+2] * 3;
 383     hw3 = iatoms[i*4+3] * 3;
 384     if (pbc == NULL)
 385     {
 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];
 392         /* 6 flops */
 393
 394         rvec_sub(after+hw2,after+ow1,doh2);
 395         rvec_sub(after+hw3,after+ow1,doh3);
 396         /* 6 flops */
 397     }
 398     else
 399     {
 400         pbc_dx_aiuc(pbc,b4+hw2,b4+ow1,dx);
 401         xb0 = dx[XX];
 402         yb0 = dx[YY];
 403         zb0 = dx[ZZ];
 404         pbc_dx_aiuc(pbc,b4+hw3,b4+ow1,dx);
 405         xc0 = dx[XX];
 406         yc0 = dx[YY];
 407         zc0 = dx[ZZ];
 408
 409         /* Tedious way of doing pbc */
 410         is = pbc_dx_aiuc(pbc,after+hw2,after+ow1,doh2);
 411         if (is == CENTRAL)
 412         {
 413             clear_rvec(sh_hw2);
 414         }
 415         else
 416         {
 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);
 421         }
 422         is = pbc_dx_aiuc(pbc,after+hw3,after+ow1,doh3);
 423         if (is == CENTRAL)
 424         {
 425             clear_rvec(sh_hw3);
 426         }
 427         else
 428         {
 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);
 433         }
 434     }
 435
 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.
 441      */
 442     xa1 = -(doh2[XX] + doh3[XX]) * wh;
 443     ya1 = -(doh2[YY] + doh3[YY]) * wh;
 444     za1 = -(doh2[ZZ] + doh3[ZZ]) * wh;
 445
 446     xcom = after[ow1 + XX] - xa1;
 447     ycom = after[ow1 + YY] - ya1;
 448     zcom = after[ow1 + ZZ] - za1;
 449
 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;
 456     /* 15 flops */
 457
 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;
 467     /* 27 flops */
 468
 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);
 472
 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;
 482     /* 24 flops */
 483
 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;
 488     /*
 489     xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
 490     ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
 491     */
 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;
 499     /* 65 flops */
 500
 501     sinphi = za1d * gmx_invsqrt(ra*ra);
 502     tmp    = 1.0 - sinphi * sinphi;
 503     if (tmp <= 0) {
 504       bOK = FALSE;
 505     } else {
 506       tmp2   = gmx_invsqrt(tmp);
 507       cosphi = tmp*tmp2;
 508       sinpsi = (zb1d - zc1d) * irc2 * tmp2;
 509       tmp2   = 1.0 - sinpsi * sinpsi;
 510       if (tmp2 <= 0) {
 511         bOK = FALSE;
 512       } else {
 513         cospsi = tmp2*gmx_invsqrt(tmp2);
 514       }
 515     }
 516     /* 46 flops */
 517
 518     if (bOK) {
 519       ya2d =  ra * cosphi;
 520       xb2d = -rc * cospsi;
 521       t1   = -rb * cosphi;
 522       t2   =  rc * sinpsi * sinphi;
 523       yb2d =  t1 - t2;
 524       yc2d =  t1 + t2;
 525       /* 7 flops */
 526
 527       /*     --- Step3  al,be,ga                      --- */
 528       alpa   = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
 529       beta   = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
 530       gama   = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
 531       al2be2 = alpa * alpa + beta * beta;
 532       tmp2   = (al2be2 - gama * gama);
 533       sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) * gmx_invsqrt(al2be2*al2be2);
 534       /* 47 flops */
 535
 536       /*  --- Step4  A3' --- */
 537       tmp2  = 1.0 - sinthe * sinthe;
 538       costhe = tmp2*gmx_invsqrt(tmp2);
 539       xa3d = -ya2d * sinthe;
 540       ya3d = ya2d * costhe;
 541       za3d = za1d;
 542       xb3d = xb2d * costhe - yb2d * sinthe;
 543       yb3d = xb2d * sinthe + yb2d * costhe;
 544       zb3d = zb1d;
 545       xc3d = -xb2d * costhe - yc2d * sinthe;
 546       yc3d = -xb2d * sinthe + yc2d * costhe;
 547       zc3d = zc1d;
 548       /* 26 flops */
 549
 550       /*    --- Step5  A3 --- */
 551       xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
 552       ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
 553       za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
 554       xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
 555       yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
 556       zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
 557       xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
 558       yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
 559       zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
 560       /* 45 flops */
 561       after[ow1] = xcom + xa3;
 562       after[ow1 + 1] = ycom + ya3;
 563       after[ow1 + 2] = zcom + za3;
 564       after[hw2] = xcom + xb3;
 565       after[hw2 + 1] = ycom + yb3;
 566       after[hw2 + 2] = zcom + zb3;
 567       after[hw3] = xcom + xc3;
 568       after[hw3 + 1] = ycom + yc3;
 569       after[hw3 + 2] = zcom + zc3;
 570       /* 9 flops */
 571
 572       if (pbc != NULL)
 573       {
 574           rvec_inc(after+hw2,sh_hw2);
 575           rvec_inc(after+hw3,sh_hw3);
 576       }
 577
 578       dax = xa3 - xa1;
 579       day = ya3 - ya1;
 580       daz = za3 - za1;
 581       dbx = xb3 - xb1;
 582       dby = yb3 - yb1;
 583       dbz = zb3 - zb1;
 584       dcx = xc3 - xc1;
 585       dcy = yc3 - yc1;
 586       dcz = zc3 - zc1;
 587       /* 9 flops, counted with the virial */
 588
 589       if (v != NULL) {
 590           v[ow1]     += dax*invdts;
 591           v[ow1 + 1] += day*invdts;
 592           v[ow1 + 2] += daz*invdts;
 593           v[hw2]     += dbx*invdts;
 594           v[hw2 + 1] += dby*invdts;
 595           v[hw2 + 2] += dbz*invdts;
 596           v[hw3]     += dcx*invdts;
 597           v[hw3 + 1] += dcy*invdts;
 598           v[hw3 + 2] += dcz*invdts;
 599           /* 3*6 flops */
 600       }
 601
 602       if (ow1 < CalcVirAtomEnd) {
 603           mdax = mOs*dax;
 604           mday = mOs*day;
 605           mdaz = mOs*daz;
 606           mdbx = mHs*dbx;
 607           mdby = mHs*dby;
 608           mdbz = mHs*dbz;
 609           mdcx = mHs*dcx;
 610           mdcy = mHs*dcy;
 611           mdcz = mHs*dcz;
 612           vir_r_m_dr[XX][XX] -= b4[ow1  ]*mdax + (b4[ow1  ]+xb0)*mdbx + (b4[ow1  ]+xc0)*mdcx;
 613           vir_r_m_dr[XX][YY] -= b4[ow1  ]*mday + (b4[ow1  ]+xb0)*mdby + (b4[ow1  ]+xc0)*mdcy;
 614           vir_r_m_dr[XX][ZZ] -= b4[ow1  ]*mdaz + (b4[ow1  ]+xb0)*mdbz + (b4[ow1  ]+xc0)*mdcz;
 615           vir_r_m_dr[YY][XX] -= b4[ow1+1]*mdax + (b4[ow1+1]+yb0)*mdbx + (b4[ow1+1]+yc0)*mdcx;
 616           vir_r_m_dr[YY][YY] -= b4[ow1+1]*mday + (b4[ow1+1]+yb0)*mdby + (b4[ow1+1]+yc0)*mdcy;
 617           vir_r_m_dr[YY][ZZ] -= b4[ow1+1]*mdaz + (b4[ow1+1]+yb0)*mdbz + (b4[ow1+1]+yc0)*mdcz;
 618           vir_r_m_dr[ZZ][XX] -= b4[ow1+2]*mdax + (b4[ow1+2]+zb0)*mdbx + (b4[ow1+2]+zc0)*mdcx;
 619           vir_r_m_dr[ZZ][YY] -= b4[ow1+2]*mday + (b4[ow1+2]+zb0)*mdby + (b4[ow1+2]+zc0)*mdcy;
 620           vir_r_m_dr[ZZ][ZZ] -= b4[ow1+2]*mdaz + (b4[ow1+2]+zb0)*mdbz + (b4[ow1+2]+zc0)*mdcz;
 621         /* 3*24 - 9 flops */
 622       }
 623     } else {
 624       *error = i;
 625     }
 626 #ifdef DEBUG
 627     if (debug)
 628     {
 629         check_cons(debug,"settle",after,ow1,hw2,hw3);
 630     }
 631 #endif
 632   }
 633 }