src/gromacs/mdlib/csettle.c

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