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 "config.h"
  40
  41 #include <math.h>
  42 #include <stdio.h>
  43
  44 #include "gromacs/legacyheaders/constr.h"
  45
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/pbcutil/ishift.h"
  48 #include "gromacs/pbcutil/pbc.h"
  49 #include "gromacs/utility/fatalerror.h"
  50 #include "gromacs/utility/smalloc.h"
  51
  52 typedef struct
  53 {
  54     real   mO;
  55     real   mH;
  56     real   wh;
  57     real   dOH;
  58     real   dHH;
  59     real   ra;
  60     real   rb;
  61     real   rc;
  62     real   irc2;
  63     /* For projection */
  64     real   imO;
  65     real   imH;
  66     real   invdOH;
  67     real   invdHH;
  68     matrix invmat;
  69 } settleparam_t;
  70
  71 typedef struct gmx_settledata
  72 {
  73     settleparam_t massw;
  74     settleparam_t mass1;
  75 } t_gmx_settledata;
  76
  77
  78 static void init_proj_matrix(settleparam_t *p,
  79                              real invmO, real invmH, real dOH, real dHH)
  80 {
  81     real   imOn, imHn;
  82     matrix mat;
  83
  84     p->imO = invmO;
  85     p->imH = invmH;
  86     /* We normalize the inverse masses with imO for the matrix inversion.
  87      * so we can keep using masses of almost zero for frozen particles,
  88      * without running out of the float range in m_inv.
  89      */
  90     imOn = 1;
  91     imHn = p->imH/p->imO;
  92
  93     /* Construct the constraint coupling matrix */
  94     mat[0][0] = imOn + imHn;
  95     mat[0][1] = imOn*(1 - 0.5*dHH*dHH/(dOH*dOH));
  96     mat[0][2] = imHn*0.5*dHH/dOH;
  97     mat[1][1] = mat[0][0];
  98     mat[1][2] = mat[0][2];
  99     mat[2][2] = imHn + imHn;
 100     mat[1][0] = mat[0][1];
 101     mat[2][0] = mat[0][2];
 102     mat[2][1] = mat[1][2];
 103
 104     m_inv(mat, p->invmat);
 105
 106     msmul(p->invmat, 1/p->imO, p->invmat);
 107
 108     p->invdOH = 1/dOH;
 109     p->invdHH = 1/dHH;
 110 }
 111
 112 static void settleparam_init(settleparam_t *p,
 113                              real mO, real mH, real invmO, real invmH,
 114                              real dOH, real dHH)
 115 {
 116     double wohh;
 117
 118     p->mO     = mO;
 119     p->mH     = mH;
 120     wohh      = mO + 2.0*mH;
 121     p->wh     = mH/wohh;
 122     p->dOH    = dOH;
 123     p->dHH    = dHH;
 124     p->rc     = dHH/2.0;
 125     p->ra     = 2.0*mH*sqrt(dOH*dOH - p->rc*p->rc)/wohh;
 126     p->rb     = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
 127     p->irc2   = 1.0/dHH;
 128
 129     /* For projection: connection matrix inversion */
 130     init_proj_matrix(p, invmO, invmH, dOH, dHH);
 131
 132     if (debug)
 133     {
 134         fprintf(debug, "wh =%g, rc = %g, ra = %g\n",
 135                 p->wh, p->rc, p->ra);
 136         fprintf(debug, "rb = %g, irc2 = %g, dHH = %g, dOH = %g\n",
 137                 p->rb, p->irc2, p->dHH, p->dOH);
 138     }
 139 }
 140
 141 gmx_settledata_t settle_init(real mO, real mH, real invmO, real invmH,
 142                              real dOH, real dHH)
 143 {
 144     gmx_settledata_t settled;
 145
 146     snew(settled, 1);
 147
 148     settleparam_init(&settled->massw, mO, mH, invmO, invmH, dOH, dHH);
 149
 150     settleparam_init(&settled->mass1, 1.0, 1.0, 1.0, 1.0, dOH, dHH);
 151
 152     return settled;
 153 }
 154
 155 #ifdef DEBUG
 156 static void check_cons(FILE *fp, char *title, real x[], int OW1, int HW2, int HW3)
 157 {
 158     rvec dOH1, dOH2, dHH;
 159     int  m;
 160
 161     for (m = 0; (m < DIM); m++)
 162     {
 163         dOH1[m] = x[OW1+m]-x[HW2+m];
 164         dOH2[m] = x[OW1+m]-x[HW3+m];
 165         dHH[m]  = x[HW2+m]-x[HW3+m];
 166     }
 167     fprintf(fp, "%10s, OW1=%3d, HW2=%3d, HW3=%3d,  dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
 168             title, OW1/DIM, HW2/DIM, HW3/DIM, norm(dOH1), norm(dOH2), norm(dHH));
 169 }
 170 #endif
 171
 172
 173 void settle_proj(gmx_settledata_t settled, int econq,
 174                  int nsettle, t_iatom iatoms[],
 175                  const t_pbc *pbc,
 176                  rvec x[],
 177                  rvec *der, rvec *derp,
 178                  int calcvir_atom_end, tensor vir_r_m_dder,
 179                  t_vetavars *vetavar)
 180 {
 181     /* Settle for projection out constraint components
 182      * of derivatives of the coordinates.
 183      * Berk Hess 2008-1-10
 184      */
 185
 186     settleparam_t *p;
 187     real           imO, imH, dOH, dHH, invdOH, invdHH;
 188     matrix         invmat;
 189     int            i, m, m2, ow1, hw2, hw3;
 190     rvec           roh2, roh3, rhh, dc, fc, fcv;
 191     rvec           derm[3], derpm[3];
 192     real           invvscale, vscale_nhc, veta;
 193     real           kfacOH, kfacHH;
 194
 195     calcvir_atom_end *= DIM;
 196
 197     if (econq == econqForce)
 198     {
 199         p = &settled->mass1;
 200     }
 201     else
 202     {
 203         p = &settled->massw;
 204     }
 205     imO    = p->imO;
 206     imH    = p->imH;
 207     copy_mat(p->invmat, invmat);
 208     dOH    = p->dOH;
 209     dHH    = p->dHH;
 210     invdOH = p->invdOH;
 211     invdHH = p->invdHH;
 212
 213     veta       = vetavar->veta;
 214     vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
 215
 216 #ifdef PRAGMAS
 217 #pragma ivdep
 218 #endif
 219
 220     for (i = 0; i < nsettle; i++)
 221     {
 222         ow1 = iatoms[i*4+1];
 223         hw2 = iatoms[i*4+2];
 224         hw3 = iatoms[i*4+3];
 225
 226
 227         for (m = 0; m < DIM; m++)
 228         {
 229             /* in the velocity case, these are the velocities, so we
 230                need to modify with the pressure control velocities! */
 231
 232             derm[0][m]  = vscale_nhc*der[ow1][m] + veta*x[ow1][m];
 233             derm[1][m]  = vscale_nhc*der[hw2][m] + veta*x[hw2][m];
 234             derm[2][m]  = vscale_nhc*der[hw3][m] + veta*x[hw3][m];
 235
 236         }
 237         /* 27 flops */
 238
 239         if (pbc == NULL)
 240         {
 241             rvec_sub(x[ow1], x[hw2], roh2);
 242             rvec_sub(x[ow1], x[hw3], roh3);
 243             rvec_sub(x[hw2], x[hw3], rhh);
 244         }
 245         else
 246         {
 247             pbc_dx_aiuc(pbc, x[ow1], x[hw2], roh2);
 248             pbc_dx_aiuc(pbc, x[ow1], x[hw3], roh3);
 249             pbc_dx_aiuc(pbc, x[hw2], x[hw3], rhh);
 250         }
 251         svmul(invdOH, roh2, roh2);
 252         svmul(invdOH, roh3, roh3);
 253         svmul(invdHH, rhh, rhh);
 254         /* 18 flops */
 255
 256         /* Determine the projections of der(modified) on the bonds */
 257         clear_rvec(dc);
 258         for (m = 0; m < DIM; m++)
 259         {
 260             dc[0] += (derm[0][m] - derm[1][m])*roh2[m];
 261             dc[1] += (derm[0][m] - derm[2][m])*roh3[m];
 262             dc[2] += (derm[1][m] - derm[2][m])*rhh [m];
 263         }
 264         /* 27 flops */
 265
 266         /* Determine the correction for the three bonds */
 267         mvmul(invmat, dc, fc);
 268         /* 15 flops */
 269
 270         /* divide velocity by vscale_nhc for determining constrained velocities, since they
 271            have not yet been multiplied */
 272         svmul(1.0/vscale_nhc, fc, fcv);
 273         /* 7? flops */
 274
 275         /* Subtract the corrections from derp */
 276         for (m = 0; m < DIM; m++)
 277         {
 278             derp[ow1][m] -= imO*( fcv[0]*roh2[m] + fcv[1]*roh3[m]);
 279             derp[hw2][m] -= imH*(-fcv[0]*roh2[m] + fcv[2]*rhh [m]);
 280             derp[hw3][m] -= imH*(-fcv[1]*roh3[m] - fcv[2]*rhh [m]);
 281         }
 282
 283         /* 45 flops */
 284
 285         if (ow1 < calcvir_atom_end)
 286         {
 287             /* Determining r \dot m der is easy,
 288              * since fc contains the mass weighted corrections for der.
 289              */
 290
 291             for (m = 0; m < DIM; m++)
 292             {
 293                 for (m2 = 0; m2 < DIM; m2++)
 294                 {
 295                     vir_r_m_dder[m][m2] +=
 296                         dOH*roh2[m]*roh2[m2]*fcv[0] +
 297                         dOH*roh3[m]*roh3[m2]*fcv[1] +
 298                         dHH*rhh [m]*rhh [m2]*fcv[2];
 299                 }
 300             }
 301         }
 302     }
 303
 304     if (calcvir_atom_end > 0)
 305     {
 306         /* Correct r_m_dder, which will be used to calcualate the virial;
 307          * we need to use the unscaled multipliers in the virial.
 308          */
 309         msmul(vir_r_m_dder, 1.0/vetavar->vscale, vir_r_m_dder);
 310     }
 311 }
 312
 313
 314 void csettle(gmx_settledata_t settled,
 315              int nsettle, t_iatom iatoms[],
 316              const t_pbc *pbc,
 317              real b4[], real after[],
 318              real invdt, real *v, int CalcVirAtomEnd,
 319              tensor vir_r_m_dr,
 320              int *error,
 321              t_vetavars *vetavar)
 322 {
 323     /* ***************************************************************** */
 324     /*                                                               ** */
 325     /*    Subroutine : setlep - reset positions of TIP3P waters      ** */
 326     /*    Author : Shuichi Miyamoto                                  ** */
 327     /*    Date of last update : Oct. 1, 1992                         ** */
 328     /*                                                               ** */
 329     /*    Reference for the SETTLE algorithm                         ** */
 330     /*           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
 331     /*                                                               ** */
 332     /* ***************************************************************** */
 333
 334     /* Initialized data */
 335     settleparam_t *p;
 336     real           wh, ra, rb, rc, irc2;
 337     real           mOs, mHs, invdts;
 338
 339     /* Local variables */
 340     real gama, beta, alpa, xcom, ycom, zcom, al2be2, tmp, tmp2;
 341     real axlng, aylng, azlng, trns11, trns21, trns31, trns12, trns22,
 342          trns32, trns13, trns23, trns33, cosphi, costhe, sinphi, sinthe,
 343          cospsi, xaksxd, yaksxd, xakszd, yakszd, zakszd, zaksxd, xaksyd,
 344          xb0, yb0, zb0, xc0, yc0, zc0, xa1;
 345     real ya1, za1, xb1, yb1;
 346     real zb1, xc1, yc1, zc1, yaksyd, zaksyd, sinpsi, xa3, ya3, za3,
 347          xb3, yb3, zb3, xc3, yc3, zc3, xb0d, yb0d, xc0d, yc0d,
 348          za1d, xb1d, yb1d, zb1d, xc1d, yc1d, zc1d, ya2d, xb2d, yb2d, yc2d,
 349          xa3d, ya3d, za3d, xb3d, yb3d, zb3d, xc3d, yc3d, zc3d;
 350     real     t1, t2;
 351     real     dax, day, daz, dbx, dby, dbz, dcx, dcy, dcz;
 352     real     mdax, mday, mdaz, mdbx, mdby, mdbz, mdcx, mdcy, mdcz;
 353
 354     gmx_bool bOK;
 355
 356     int      i, ow1, hw2, hw3;
 357
 358     rvec     dx, sh_hw2 = {0, 0, 0}, sh_hw3 = {0, 0, 0};
 359     rvec     doh2, doh3;
 360     int      is;
 361
 362     *error = -1;
 363
 364     CalcVirAtomEnd *= 3;
 365
 366     p    = &settled->massw;
 367     wh   = p->wh;
 368     rc   = p->rc;
 369     ra   = p->ra;
 370     rb   = p->rb;
 371     irc2 = p->irc2;
 372
 373     mOs    = p->mO / vetavar->rvscale;
 374     mHs    = p->mH / vetavar->rvscale;
 375     invdts = invdt / vetavar->rscale;
 376
 377 #ifdef PRAGMAS
 378 #pragma ivdep
 379 #endif
 380     for (i = 0; i < nsettle; ++i)
 381     {
 382         bOK = TRUE;
 383         /*    --- Step1  A1' ---      */
 384         ow1 = iatoms[i*4+1] * 3;
 385         hw2 = iatoms[i*4+2] * 3;
 386         hw3 = iatoms[i*4+3] * 3;
 387         if (pbc == NULL)
 388         {
 389             xb0 = b4[hw2 + XX] - b4[ow1 + XX];
 390             yb0 = b4[hw2 + YY] - b4[ow1 + YY];
 391             zb0 = b4[hw2 + ZZ] - b4[ow1 + ZZ];
 392             xc0 = b4[hw3 + XX] - b4[ow1 + XX];
 393             yc0 = b4[hw3 + YY] - b4[ow1 + YY];
 394             zc0 = b4[hw3 + ZZ] - b4[ow1 + ZZ];
 395             /* 6 flops */
 396
 397             rvec_sub(after+hw2, after+ow1, doh2);
 398             rvec_sub(after+hw3, after+ow1, doh3);
 399             /* 6 flops */
 400         }
 401         else
 402         {
 403             pbc_dx_aiuc(pbc, b4+hw2, b4+ow1, dx);
 404             xb0 = dx[XX];
 405             yb0 = dx[YY];
 406             zb0 = dx[ZZ];
 407             pbc_dx_aiuc(pbc, b4+hw3, b4+ow1, dx);
 408             xc0 = dx[XX];
 409             yc0 = dx[YY];
 410             zc0 = dx[ZZ];
 411
 412             /* Tedious way of doing pbc */
 413             is = pbc_dx_aiuc(pbc, after+hw2, after+ow1, doh2);
 414             if (is == CENTRAL)
 415             {
 416                 clear_rvec(sh_hw2);
 417             }
 418             else
 419             {
 420                 sh_hw2[XX] = after[hw2 + XX] - (after[ow1 + XX] + doh2[XX]);
 421                 sh_hw2[YY] = after[hw2 + YY] - (after[ow1 + YY] + doh2[YY]);
 422                 sh_hw2[ZZ] = after[hw2 + ZZ] - (after[ow1 + ZZ] + doh2[ZZ]);
 423                 rvec_dec(after+hw2, sh_hw2);
 424             }
 425             is = pbc_dx_aiuc(pbc, after+hw3, after+ow1, doh3);
 426             if (is == CENTRAL)
 427             {
 428                 clear_rvec(sh_hw3);
 429             }
 430             else
 431             {
 432                 sh_hw3[XX] = after[hw3 + XX] - (after[ow1 + XX] + doh3[XX]);
 433                 sh_hw3[YY] = after[hw3 + YY] - (after[ow1 + YY] + doh3[YY]);
 434                 sh_hw3[ZZ] = after[hw3 + ZZ] - (after[ow1 + ZZ] + doh3[ZZ]);
 435                 rvec_dec(after+hw3, sh_hw3);
 436             }
 437         }
 438
 439         /* Not calculating the center of mass using the oxygen position
 440          * and the O-H distances, as done below, will make SETTLE
 441          * the largest source of energy drift for simulations of water,
 442          * as then the oxygen coordinate is multiplied by 0.89 at every step,
 443          * which can then transfer a systematic rounding to the oxygen velocity.
 444          */
 445         xa1 = -(doh2[XX] + doh3[XX]) * wh;
 446         ya1 = -(doh2[YY] + doh3[YY]) * wh;
 447         za1 = -(doh2[ZZ] + doh3[ZZ]) * wh;
 448
 449         xcom = after[ow1 + XX] - xa1;
 450         ycom = after[ow1 + YY] - ya1;
 451         zcom = after[ow1 + ZZ] - za1;
 452
 453         xb1 = after[hw2 + XX] - xcom;
 454         yb1 = after[hw2 + YY] - ycom;
 455         zb1 = after[hw2 + ZZ] - zcom;
 456         xc1 = after[hw3 + XX] - xcom;
 457         yc1 = after[hw3 + YY] - ycom;
 458         zc1 = after[hw3 + ZZ] - zcom;
 459         /* 15 flops */
 460
 461         xakszd = yb0 * zc0 - zb0 * yc0;
 462         yakszd = zb0 * xc0 - xb0 * zc0;
 463         zakszd = xb0 * yc0 - yb0 * xc0;
 464         xaksxd = ya1 * zakszd - za1 * yakszd;
 465         yaksxd = za1 * xakszd - xa1 * zakszd;
 466         zaksxd = xa1 * yakszd - ya1 * xakszd;
 467         xaksyd = yakszd * zaksxd - zakszd * yaksxd;
 468         yaksyd = zakszd * xaksxd - xakszd * zaksxd;
 469         zaksyd = xakszd * yaksxd - yakszd * xaksxd;
 470         /* 27 flops */
 471
 472         axlng = gmx_invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
 473         aylng = gmx_invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
 474         azlng = gmx_invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
 475
 476         trns11 = xaksxd * axlng;
 477         trns21 = yaksxd * axlng;
 478         trns31 = zaksxd * axlng;
 479         trns12 = xaksyd * aylng;
 480         trns22 = yaksyd * aylng;
 481         trns32 = zaksyd * aylng;
 482         trns13 = xakszd * azlng;
 483         trns23 = yakszd * azlng;
 484         trns33 = zakszd * azlng;
 485         /* 24 flops */
 486
 487         xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
 488         yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
 489         xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
 490         yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
 491         /*
 492            xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
 493            ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
 494          */
 495         za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
 496         xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
 497         yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
 498         zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
 499         xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
 500         yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
 501         zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
 502         /* 65 flops */
 503
 504         sinphi = za1d * gmx_invsqrt(ra*ra);
 505         tmp    = 1.0 - sinphi * sinphi;
 506         if (tmp <= 0)
 507         {
 508             bOK = FALSE;
 509         }
 510         else
 511         {
 512             tmp2   = gmx_invsqrt(tmp);
 513             cosphi = tmp*tmp2;
 514             sinpsi = (zb1d - zc1d) * irc2 * tmp2;
 515             tmp2   = 1.0 - sinpsi * sinpsi;
 516             if (tmp2 <= 0)
 517             {
 518                 bOK = FALSE;
 519             }
 520             else
 521             {
 522                 cospsi = tmp2*gmx_invsqrt(tmp2);
 523             }
 524         }
 525         /* 46 flops */
 526
 527         if (bOK)
 528         {
 529             ya2d =  ra * cosphi;
 530             xb2d = -rc * cospsi;
 531             t1   = -rb * cosphi;
 532             t2   =  rc * sinpsi * sinphi;
 533             yb2d =  t1 - t2;
 534             yc2d =  t1 + t2;
 535             /* 7 flops */
 536
 537             /*     --- Step3  al,be,ga            --- */
 538             alpa   = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
 539             beta   = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
 540             gama   = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
 541             al2be2 = alpa * alpa + beta * beta;
 542             tmp2   = (al2be2 - gama * gama);
 543             sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) * gmx_invsqrt(al2be2*al2be2);
 544             /* 47 flops */
 545
 546             /*  --- Step4  A3' --- */
 547             tmp2   = 1.0 - sinthe * sinthe;
 548             costhe = tmp2*gmx_invsqrt(tmp2);
 549             xa3d   = -ya2d * sinthe;
 550             ya3d   = ya2d * costhe;
 551             za3d   = za1d;
 552             xb3d   = xb2d * costhe - yb2d * sinthe;
 553             yb3d   = xb2d * sinthe + yb2d * costhe;
 554             zb3d   = zb1d;
 555             xc3d   = -xb2d * costhe - yc2d * sinthe;
 556             yc3d   = -xb2d * sinthe + yc2d * costhe;
 557             zc3d   = zc1d;
 558             /* 26 flops */
 559
 560             /*    --- Step5  A3 --- */
 561             xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
 562             ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
 563             za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
 564             xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
 565             yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
 566             zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
 567             xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
 568             yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
 569             zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
 570             /* 45 flops */
 571             after[ow1]     = xcom + xa3;
 572             after[ow1 + 1] = ycom + ya3;
 573             after[ow1 + 2] = zcom + za3;
 574             after[hw2]     = xcom + xb3;
 575             after[hw2 + 1] = ycom + yb3;
 576             after[hw2 + 2] = zcom + zb3;
 577             after[hw3]     = xcom + xc3;
 578             after[hw3 + 1] = ycom + yc3;
 579             after[hw3 + 2] = zcom + zc3;
 580             /* 9 flops */
 581
 582             if (pbc != NULL)
 583             {
 584                 rvec_inc(after+hw2, sh_hw2);
 585                 rvec_inc(after+hw3, sh_hw3);
 586             }
 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 != NULL)
 600             {
 601                 v[ow1]     += dax*invdts;
 602                 v[ow1 + 1] += day*invdts;
 603                 v[ow1 + 2] += daz*invdts;
 604                 v[hw2]     += dbx*invdts;
 605                 v[hw2 + 1] += dby*invdts;
 606                 v[hw2 + 2] += dbz*invdts;
 607                 v[hw3]     += dcx*invdts;
 608                 v[hw3 + 1] += dcy*invdts;
 609                 v[hw3 + 2] += dcz*invdts;
 610                 /* 3*6 flops */
 611             }
 612
 613             if (ow1 < CalcVirAtomEnd)
 614             {
 615                 mdax                = mOs*dax;
 616                 mday                = mOs*day;
 617                 mdaz                = mOs*daz;
 618                 mdbx                = mHs*dbx;
 619                 mdby                = mHs*dby;
 620                 mdbz                = mHs*dbz;
 621                 mdcx                = mHs*dcx;
 622                 mdcy                = mHs*dcy;
 623                 mdcz                = mHs*dcz;
 624                 vir_r_m_dr[XX][XX] -= b4[ow1  ]*mdax + (b4[ow1  ]+xb0)*mdbx + (b4[ow1  ]+xc0)*mdcx;
 625                 vir_r_m_dr[XX][YY] -= b4[ow1  ]*mday + (b4[ow1  ]+xb0)*mdby + (b4[ow1  ]+xc0)*mdcy;
 626                 vir_r_m_dr[XX][ZZ] -= b4[ow1  ]*mdaz + (b4[ow1  ]+xb0)*mdbz + (b4[ow1  ]+xc0)*mdcz;
 627                 vir_r_m_dr[YY][XX] -= b4[ow1+1]*mdax + (b4[ow1+1]+yb0)*mdbx + (b4[ow1+1]+yc0)*mdcx;
 628                 vir_r_m_dr[YY][YY] -= b4[ow1+1]*mday + (b4[ow1+1]+yb0)*mdby + (b4[ow1+1]+yc0)*mdcy;
 629                 vir_r_m_dr[YY][ZZ] -= b4[ow1+1]*mdaz + (b4[ow1+1]+yb0)*mdbz + (b4[ow1+1]+yc0)*mdcz;
 630                 vir_r_m_dr[ZZ][XX] -= b4[ow1+2]*mdax + (b4[ow1+2]+zb0)*mdbx + (b4[ow1+2]+zc0)*mdcx;
 631                 vir_r_m_dr[ZZ][YY] -= b4[ow1+2]*mday + (b4[ow1+2]+zb0)*mdby + (b4[ow1+2]+zc0)*mdcy;
 632                 vir_r_m_dr[ZZ][ZZ] -= b4[ow1+2]*mdaz + (b4[ow1+2]+zb0)*mdbz + (b4[ow1+2]+zc0)*mdcz;
 633                 /* 3*24 - 9 flops */
 634             }
 635         }
 636         else
 637         {
 638             *error = i;
 639         }
 640 #ifdef DEBUG
 641         if (debug)
 642         {
 643             check_cons(debug, "settle", after, ow1, hw2, hw3);
 644         }
 645 #endif
 646     }
 647 }