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
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/pbcutil/ishift.h"
  46 #include "gromacs/pbcutil/pbc.h"
  47 #include "gromacs/utility/fatalerror.h"
  48 #include "gromacs/utility/smalloc.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     {
 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(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     {
 380         bOK = TRUE;
 381         /*    --- Step1  A1' ---      */
 382         ow1 = iatoms[i*4+1] * 3;
 383         hw2 = iatoms[i*4+2] * 3;
 384         hw3 = iatoms[i*4+3] * 3;
 385         if (pbc == NULL)
 386         {
 387             xb0 = b4[hw2 + XX] - b4[ow1 + XX];
 388             yb0 = b4[hw2 + YY] - b4[ow1 + YY];
 389             zb0 = b4[hw2 + ZZ] - b4[ow1 + ZZ];
 390             xc0 = b4[hw3 + XX] - b4[ow1 + XX];
 391             yc0 = b4[hw3 + YY] - b4[ow1 + YY];
 392             zc0 = b4[hw3 + ZZ] - b4[ow1 + ZZ];
 393             /* 6 flops */
 394
 395             rvec_sub(after+hw2, after+ow1, doh2);
 396             rvec_sub(after+hw3, after+ow1, doh3);
 397             /* 6 flops */
 398         }
 399         else
 400         {
 401             pbc_dx_aiuc(pbc, b4+hw2, b4+ow1, dx);
 402             xb0 = dx[XX];
 403             yb0 = dx[YY];
 404             zb0 = dx[ZZ];
 405             pbc_dx_aiuc(pbc, b4+hw3, b4+ow1, dx);
 406             xc0 = dx[XX];
 407             yc0 = dx[YY];
 408             zc0 = dx[ZZ];
 409
 410             /* Tedious way of doing pbc */
 411             is = pbc_dx_aiuc(pbc, after+hw2, after+ow1, doh2);
 412             if (is == CENTRAL)
 413             {
 414                 clear_rvec(sh_hw2);
 415             }
 416             else
 417             {
 418                 sh_hw2[XX] = after[hw2 + XX] - (after[ow1 + XX] + doh2[XX]);
 419                 sh_hw2[YY] = after[hw2 + YY] - (after[ow1 + YY] + doh2[YY]);
 420                 sh_hw2[ZZ] = after[hw2 + ZZ] - (after[ow1 + ZZ] + doh2[ZZ]);
 421                 rvec_dec(after+hw2, sh_hw2);
 422             }
 423             is = pbc_dx_aiuc(pbc, after+hw3, after+ow1, doh3);
 424             if (is == CENTRAL)
 425             {
 426                 clear_rvec(sh_hw3);
 427             }
 428             else
 429             {
 430                 sh_hw3[XX] = after[hw3 + XX] - (after[ow1 + XX] + doh3[XX]);
 431                 sh_hw3[YY] = after[hw3 + YY] - (after[ow1 + YY] + doh3[YY]);
 432                 sh_hw3[ZZ] = after[hw3 + ZZ] - (after[ow1 + ZZ] + doh3[ZZ]);
 433                 rvec_dec(after+hw3, sh_hw3);
 434             }
 435         }
 436
 437         /* Not calculating the center of mass using the oxygen position
 438          * and the O-H distances, as done below, will make SETTLE
 439          * the largest source of energy drift for simulations of water,
 440          * as then the oxygen coordinate is multiplied by 0.89 at every step,
 441          * which can then transfer a systematic rounding to the oxygen velocity.
 442          */
 443         xa1 = -(doh2[XX] + doh3[XX]) * wh;
 444         ya1 = -(doh2[YY] + doh3[YY]) * wh;
 445         za1 = -(doh2[ZZ] + doh3[ZZ]) * wh;
 446
 447         xcom = after[ow1 + XX] - xa1;
 448         ycom = after[ow1 + YY] - ya1;
 449         zcom = after[ow1 + ZZ] - za1;
 450
 451         xb1 = after[hw2 + XX] - xcom;
 452         yb1 = after[hw2 + YY] - ycom;
 453         zb1 = after[hw2 + ZZ] - zcom;
 454         xc1 = after[hw3 + XX] - xcom;
 455         yc1 = after[hw3 + YY] - ycom;
 456         zc1 = after[hw3 + ZZ] - zcom;
 457         /* 15 flops */
 458
 459         xakszd = yb0 * zc0 - zb0 * yc0;
 460         yakszd = zb0 * xc0 - xb0 * zc0;
 461         zakszd = xb0 * yc0 - yb0 * xc0;
 462         xaksxd = ya1 * zakszd - za1 * yakszd;
 463         yaksxd = za1 * xakszd - xa1 * zakszd;
 464         zaksxd = xa1 * yakszd - ya1 * xakszd;
 465         xaksyd = yakszd * zaksxd - zakszd * yaksxd;
 466         yaksyd = zakszd * xaksxd - xakszd * zaksxd;
 467         zaksyd = xakszd * yaksxd - yakszd * xaksxd;
 468         /* 27 flops */
 469
 470         axlng = gmx_invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
 471         aylng = gmx_invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
 472         azlng = gmx_invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
 473
 474         trns11 = xaksxd * axlng;
 475         trns21 = yaksxd * axlng;
 476         trns31 = zaksxd * axlng;
 477         trns12 = xaksyd * aylng;
 478         trns22 = yaksyd * aylng;
 479         trns32 = zaksyd * aylng;
 480         trns13 = xakszd * azlng;
 481         trns23 = yakszd * azlng;
 482         trns33 = zakszd * azlng;
 483         /* 24 flops */
 484
 485         xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
 486         yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
 487         xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
 488         yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
 489         /*
 490            xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
 491            ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
 492          */
 493         za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
 494         xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
 495         yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
 496         zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
 497         xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
 498         yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
 499         zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
 500         /* 65 flops */
 501
 502         sinphi = za1d * gmx_invsqrt(ra*ra);
 503         tmp    = 1.0 - sinphi * sinphi;
 504         if (tmp <= 0)
 505         {
 506             bOK = FALSE;
 507         }
 508         else
 509         {
 510             tmp2   = gmx_invsqrt(tmp);
 511             cosphi = tmp*tmp2;
 512             sinpsi = (zb1d - zc1d) * irc2 * tmp2;
 513             tmp2   = 1.0 - sinpsi * sinpsi;
 514             if (tmp2 <= 0)
 515             {
 516                 bOK = FALSE;
 517             }
 518             else
 519             {
 520                 cospsi = tmp2*gmx_invsqrt(tmp2);
 521             }
 522         }
 523         /* 46 flops */
 524
 525         if (bOK)
 526         {
 527             ya2d =  ra * cosphi;
 528             xb2d = -rc * cospsi;
 529             t1   = -rb * cosphi;
 530             t2   =  rc * sinpsi * sinphi;
 531             yb2d =  t1 - t2;
 532             yc2d =  t1 + t2;
 533             /* 7 flops */
 534
 535             /*     --- Step3  al,be,ga            --- */
 536             alpa   = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
 537             beta   = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
 538             gama   = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
 539             al2be2 = alpa * alpa + beta * beta;
 540             tmp2   = (al2be2 - gama * gama);
 541             sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) * gmx_invsqrt(al2be2*al2be2);
 542             /* 47 flops */
 543
 544             /*  --- Step4  A3' --- */
 545             tmp2   = 1.0 - sinthe * sinthe;
 546             costhe = tmp2*gmx_invsqrt(tmp2);
 547             xa3d   = -ya2d * sinthe;
 548             ya3d   = ya2d * costhe;
 549             za3d   = za1d;
 550             xb3d   = xb2d * costhe - yb2d * sinthe;
 551             yb3d   = xb2d * sinthe + yb2d * costhe;
 552             zb3d   = zb1d;
 553             xc3d   = -xb2d * costhe - yc2d * sinthe;
 554             yc3d   = -xb2d * sinthe + yc2d * costhe;
 555             zc3d   = zc1d;
 556             /* 26 flops */
 557
 558             /*    --- Step5  A3 --- */
 559             xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
 560             ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
 561             za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
 562             xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
 563             yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
 564             zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
 565             xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
 566             yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
 567             zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
 568             /* 45 flops */
 569             after[ow1]     = xcom + xa3;
 570             after[ow1 + 1] = ycom + ya3;
 571             after[ow1 + 2] = zcom + za3;
 572             after[hw2]     = xcom + xb3;
 573             after[hw2 + 1] = ycom + yb3;
 574             after[hw2 + 2] = zcom + zb3;
 575             after[hw3]     = xcom + xc3;
 576             after[hw3 + 1] = ycom + yc3;
 577             after[hw3 + 2] = zcom + zc3;
 578             /* 9 flops */
 579
 580             if (pbc != NULL)
 581             {
 582                 rvec_inc(after+hw2, sh_hw2);
 583                 rvec_inc(after+hw3, sh_hw3);
 584             }
 585
 586             dax = xa3 - xa1;
 587             day = ya3 - ya1;
 588             daz = za3 - za1;
 589             dbx = xb3 - xb1;
 590             dby = yb3 - yb1;
 591             dbz = zb3 - zb1;
 592             dcx = xc3 - xc1;
 593             dcy = yc3 - yc1;
 594             dcz = zc3 - zc1;
 595             /* 9 flops, counted with the virial */
 596
 597             if (v != NULL)
 598             {
 599                 v[ow1]     += dax*invdts;
 600                 v[ow1 + 1] += day*invdts;
 601                 v[ow1 + 2] += daz*invdts;
 602                 v[hw2]     += dbx*invdts;
 603                 v[hw2 + 1] += dby*invdts;
 604                 v[hw2 + 2] += dbz*invdts;
 605                 v[hw3]     += dcx*invdts;
 606                 v[hw3 + 1] += dcy*invdts;
 607                 v[hw3 + 2] += dcz*invdts;
 608                 /* 3*6 flops */
 609             }
 610
 611             if (ow1 < CalcVirAtomEnd)
 612             {
 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                 vir_r_m_dr[XX][XX] -= b4[ow1  ]*mdax + (b4[ow1  ]+xb0)*mdbx + (b4[ow1  ]+xc0)*mdcx;
 623                 vir_r_m_dr[XX][YY] -= b4[ow1  ]*mday + (b4[ow1  ]+xb0)*mdby + (b4[ow1  ]+xc0)*mdcy;
 624                 vir_r_m_dr[XX][ZZ] -= b4[ow1  ]*mdaz + (b4[ow1  ]+xb0)*mdbz + (b4[ow1  ]+xc0)*mdcz;
 625                 vir_r_m_dr[YY][XX] -= b4[ow1+1]*mdax + (b4[ow1+1]+yb0)*mdbx + (b4[ow1+1]+yc0)*mdcx;
 626                 vir_r_m_dr[YY][YY] -= b4[ow1+1]*mday + (b4[ow1+1]+yb0)*mdby + (b4[ow1+1]+yc0)*mdcy;
 627                 vir_r_m_dr[YY][ZZ] -= b4[ow1+1]*mdaz + (b4[ow1+1]+yb0)*mdbz + (b4[ow1+1]+yc0)*mdcz;
 628                 vir_r_m_dr[ZZ][XX] -= b4[ow1+2]*mdax + (b4[ow1+2]+zb0)*mdbx + (b4[ow1+2]+zc0)*mdcx;
 629                 vir_r_m_dr[ZZ][YY] -= b4[ow1+2]*mday + (b4[ow1+2]+zb0)*mdby + (b4[ow1+2]+zc0)*mdcy;
 630                 vir_r_m_dr[ZZ][ZZ] -= b4[ow1+2]*mdaz + (b4[ow1+2]+zb0)*mdbz + (b4[ow1+2]+zc0)*mdcz;
 631                 /* 3*24 - 9 flops */
 632             }
 633         }
 634         else
 635         {
 636             *error = i;
 637         }
 638 #ifdef DEBUG
 639         if (debug)
 640         {
 641             check_cons(debug, "settle", after, ow1, hw2, hw3);
 642         }
 643 #endif
 644     }
 645 }