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