src/gromacs/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  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <math.h>
  40 #include <stdio.h>
  41
  42 #include "gromacs/legacyheaders/constr.h"
  43 #include "gromacs/math/vec.h"
  44 #include "gromacs/pbcutil/ishift.h"
  45 #include "gromacs/pbcutil/pbc.h"
  46 #include "gromacs/utility/fatalerror.h"
  47 #include "gromacs/utility/smalloc.h"
  48
  49 typedef struct
  50 {
  51     real   mO;
  52     real   mH;
  53     real   wh;
  54     real   dOH;
  55     real   dHH;
  56     real   ra;
  57     real   rb;
  58     real   rc;
  59     real   irc2;
  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     double wohh;
 114
 115     p->mO     = mO;
 116     p->mH     = mH;
 117     wohh      = mO + 2.0*mH;
 118     p->wh     = mH/wohh;
 119     p->dOH    = dOH;
 120     p->dHH    = dHH;
 121     p->rc     = dHH/2.0;
 122     p->ra     = 2.0*mH*sqrt(dOH*dOH - p->rc*p->rc)/wohh;
 123     p->rb     = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
 124     p->irc2   = 1.0/dHH;
 125
 126     /* For projection: connection matrix inversion */
 127     init_proj_matrix(p, invmO, invmH, dOH, dHH);
 128
 129     if (debug)
 130     {
 131         fprintf(debug, "wh =%g, rc = %g, ra = %g\n",
 132                 p->wh, p->rc, p->ra);
 133         fprintf(debug, "rb = %g, irc2 = %g, dHH = %g, dOH = %g\n",
 134                 p->rb, p->irc2, p->dHH, p->dOH);
 135     }
 136 }
 137
 138 gmx_settledata_t settle_init(real mO, real mH, real invmO, real invmH,
 139                              real dOH, real dHH)
 140 {
 141     gmx_settledata_t settled;
 142
 143     snew(settled, 1);
 144
 145     settleparam_init(&settled->massw, mO, mH, invmO, invmH, dOH, dHH);
 146
 147     settleparam_init(&settled->mass1, 1.0, 1.0, 1.0, 1.0, dOH, dHH);
 148
 149     return settled;
 150 }
 151
 152 #ifdef DEBUG
 153 static void check_cons(FILE *fp, char *title, real x[], int OW1, int HW2, int HW3)
 154 {
 155     rvec dOH1, dOH2, dHH;
 156     int  m;
 157
 158     for (m = 0; (m < DIM); m++)
 159     {
 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(gmx_settledata_t settled, int econq,
 171                  int nsettle, t_iatom iatoms[],
 172                  const t_pbc *pbc,
 173                  rvec x[],
 174                  rvec *der, rvec *derp,
 175                  int calcvir_atom_end, tensor vir_r_m_dder,
 176                  t_vetavars *vetavar)
 177 {
 178     /* Settle for projection out constraint components
 179      * of derivatives of the coordinates.
 180      * Berk Hess 2008-1-10
 181      */
 182
 183     settleparam_t *p;
 184     real           imO, imH, dOH, dHH, invdOH, invdHH;
 185     matrix         invmat;
 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;
 190     real           kfacOH, kfacHH;
 191
 192     calcvir_atom_end *= DIM;
 193
 194     if (econq == econqForce)
 195     {
 196         p = &settled->mass1;
 197     }
 198     else
 199     {
 200         p = &settled->massw;
 201     }
 202     imO    = p->imO;
 203     imH    = p->imH;
 204     copy_mat(p->invmat, invmat);
 205     dOH    = p->dOH;
 206     dHH    = p->dHH;
 207     invdOH = p->invdOH;
 208     invdHH = p->invdHH;
 209
 210     veta       = vetavar->veta;
 211     vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
 212
 213 #ifdef PRAGMAS
 214 #pragma ivdep
 215 #endif
 216
 217     for (i = 0; i < nsettle; i++)
 218     {
 219         ow1 = iatoms[i*4+1];
 220         hw2 = iatoms[i*4+2];
 221         hw3 = iatoms[i*4+3];
 222
 223
 224         for (m = 0; m < DIM; m++)
 225         {
 226             /* in the velocity case, these are the velocities, so we
 227                need to modify with the pressure control velocities! */
 228
 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];
 232
 233         }
 234         /* 27 flops */
 235
 236         if (pbc == NULL)
 237         {
 238             rvec_sub(x[ow1], x[hw2], roh2);
 239             rvec_sub(x[ow1], x[hw3], roh3);
 240             rvec_sub(x[hw2], x[hw3], rhh);
 241         }
 242         else
 243         {
 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);
 247         }
 248         svmul(invdOH, roh2, roh2);
 249         svmul(invdOH, roh3, roh3);
 250         svmul(invdHH, rhh, rhh);
 251         /* 18 flops */
 252
 253         /* Determine the projections of der(modified) on the bonds */
 254         clear_rvec(dc);
 255         for (m = 0; m < DIM; m++)
 256         {
 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];
 260         }
 261         /* 27 flops */
 262
 263         /* Determine the correction for the three bonds */
 264         mvmul(invmat, dc, fc);
 265         /* 15 flops */
 266
 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);
 270         /* 7? flops */
 271
 272         /* Subtract the corrections from derp */
 273         for (m = 0; m < DIM; m++)
 274         {
 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]);
 278         }
 279
 280         /* 45 flops */
 281
 282         if (ow1 < calcvir_atom_end)
 283         {
 284             /* Determining r \dot m der is easy,
 285              * since fc contains the mass weighted corrections for der.
 286              */
 287
 288             for (m = 0; m < DIM; m++)
 289             {
 290                 for (m2 = 0; m2 < DIM; m2++)
 291                 {
 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];
 296                 }
 297             }
 298         }
 299     }
 300
 301     if (calcvir_atom_end > 0)
 302     {
 303         /* Correct r_m_dder, which will be used to calcualate the virial;
 304          * we need to use the unscaled multipliers in the virial.
 305          */
 306         msmul(vir_r_m_dder, 1.0/vetavar->vscale, vir_r_m_dder);
 307     }
 308 }
 309
 310
 311 void csettle(gmx_settledata_t settled,
 312              int nsettle, t_iatom iatoms[],
 313              const t_pbc *pbc,
 314              real b4[], real after[],
 315              real invdt, real *v, int CalcVirAtomEnd,
 316              tensor vir_r_m_dr,
 317              int *error,
 318              t_vetavars *vetavar)
 319 {
 320     /* ***************************************************************** */
 321     /*                                                               ** */
 322     /*    Subroutine : setlep - reset positions of TIP3P waters      ** */
 323     /*    Author : Shuichi Miyamoto                                  ** */
 324     /*    Date of last update : Oct. 1, 1992                         ** */
 325     /*                                                               ** */
 326     /*    Reference for the SETTLE algorithm                         ** */
 327     /*           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
 328     /*                                                               ** */
 329     /* ***************************************************************** */
 330
 331     /* Initialized data */
 332     settleparam_t *p;
 333     real           wh, ra, rb, rc, irc2;
 334     real           mOs, mHs, invdts;
 335
 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;
 347     real     t1, t2;
 348     real     dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
 349     real     mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
 350
 351     gmx_bool bOK;
 352
 353     int      i, ow1, hw2, hw3;
 354
 355     rvec     dx, sh_hw2 = {0, 0, 0}, sh_hw3 = {0, 0, 0};
 356     rvec     doh2, doh3;
 357     int      is;
 358
 359     *error = -1;
 360
 361     CalcVirAtomEnd *= 3;
 362
 363     p    = &settled->massw;
 364     wh   = p->wh;
 365     rc   = p->rc;
 366     ra   = p->ra;
 367     rb   = p->rb;
 368     irc2 = p->irc2;
 369
 370     mOs    = p->mO / vetavar->rvscale;
 371     mHs    = p->mH / vetavar->rvscale;
 372     invdts = invdt / vetavar->rscale;
 373
 374 #ifdef PRAGMAS
 375 #pragma ivdep
 376 #endif
 377     for (i = 0; i < nsettle; ++i)
 378     {
 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         {
 505             bOK = FALSE;
 506         }
 507         else
 508         {
 509             tmp2   = gmx_invsqrt(tmp);
 510             cosphi = tmp*tmp2;
 511             sinpsi = (zb1d - zc1d) * irc2 * tmp2;
 512             tmp2   = 1.0 - sinpsi * sinpsi;
 513             if (tmp2 <= 0)
 514             {
 515                 bOK = FALSE;
 516             }
 517             else
 518             {
 519                 cospsi = tmp2*gmx_invsqrt(tmp2);
 520             }
 521         }
 522         /* 46 flops */
 523
 524         if (bOK)
 525         {
 526             ya2d =  ra * cosphi;
 527             xb2d = -rc * cospsi;
 528             t1   = -rb * cosphi;
 529             t2   =  rc * sinpsi * sinphi;
 530             yb2d =  t1 - t2;
 531             yc2d =  t1 + t2;
 532             /* 7 flops */
 533
 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);
 541             /* 47 flops */
 542
 543             /*  --- Step4  A3' --- */
 544             tmp2   = 1.0 - sinthe * sinthe;
 545             costhe = tmp2*gmx_invsqrt(tmp2);
 546             xa3d   = -ya2d * sinthe;
 547             ya3d   = ya2d * costhe;
 548             za3d   = za1d;
 549             xb3d   = xb2d * costhe - yb2d * sinthe;
 550             yb3d   = xb2d * sinthe + yb2d * costhe;
 551             zb3d   = zb1d;
 552             xc3d   = -xb2d * costhe - yc2d * sinthe;
 553             yc3d   = -xb2d * sinthe + yc2d * costhe;
 554             zc3d   = zc1d;
 555             /* 26 flops */
 556
 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;
 567             /* 45 flops */
 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;
 577             /* 9 flops */
 578
 579             if (pbc != NULL)
 580             {
 581                 rvec_inc(after+hw2, sh_hw2);
 582                 rvec_inc(after+hw3, sh_hw3);
 583             }
 584
 585             dax = xa3 - xa1;
 586             day = ya3 - ya1;
 587             daz = za3 - za1;
 588             dbx = xb3 - xb1;
 589             dby = yb3 - yb1;
 590             dbz = zb3 - zb1;
 591             dcx = xc3 - xc1;
 592             dcy = yc3 - yc1;
 593             dcz = zc3 - zc1;
 594             /* 9 flops, counted with the virial */
 595
 596             if (v != NULL)
 597             {
 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;
 607                 /* 3*6 flops */
 608             }
 609
 610             if (ow1 < CalcVirAtomEnd)
 611             {
 612                 mdax                = mOs*dax;
 613                 mday                = mOs*day;
 614                 mdaz                = mOs*daz;
 615                 mdbx                = mHs*dbx;
 616                 mdby                = mHs*dby;
 617                 mdbz                = mHs*dbz;
 618                 mdcx                = mHs*dcx;
 619                 mdcy                = mHs*dcy;
 620                 mdcz                = mHs*dcz;
 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;
 630                 /* 3*24 - 9 flops */
 631             }
 632         }
 633         else
 634         {
 635             *error = i;
 636         }
 637 #ifdef DEBUG
 638         if (debug)
 639         {
 640             check_cons(debug, "settle", after, ow1, hw2, hw3);
 641         }
 642 #endif
 643     }
 644 }