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     double wo;
  53     double wh;
  54     double wohh;
  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     p->mO     = mO;
 115     p->mH     = mH;
 116     p->wo     = mO;
 117     p->wh     = mH;
 118     p->wohh   = mO + 2.0*mH;
 119     p->dOH    = dOH;
 120     p->dHH    = dHH;
 121     p->rc     = dHH/2.0;
 122     p->ra     = 2.0*p->wh*sqrt(dOH*dOH - p->rc*p->rc)/p->wohh;
 123     p->rb     = sqrt(dOH*dOH - p->rc*p->rc) - p->ra;
 124     p->irc2   = 1.0/dHH;
 125
 126     p->wo    /= p->wohh;
 127     p->wh    /= p->wohh;
 128
 129     /* For projection: connection matrix inversion */
 130     init_proj_matrix(p,invmO,invmH,dOH,dHH);
 131
 132     if (debug)
 133     {
 134         fprintf(debug,"wo = %g, wh =%g, wohh = %g, rc = %g, ra = %g\n",
 135                 p->wo,p->wh,p->wohh,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     dOH1[m]=x[OW1+m]-x[HW2+m];
 163     dOH2[m]=x[OW1+m]-x[HW3+m];
 164     dHH[m]=x[HW2+m]-x[HW3+m];
 165   }
 166   fprintf(fp,"%10s, OW1=%3d, HW2=%3d, HW3=%3d,  dOH1: %8.3f, dOH2: %8.3f, dHH: %8.3f\n",
 167           title,OW1/DIM,HW2/DIM,HW3/DIM,norm(dOH1),norm(dOH2),norm(dHH));
 168 }
 169 #endif
 170
 171
 172 void settle_proj(FILE *fp,
 173                  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 rmdder,t_vetavars *vetavar)
 179 {
 180     /* Settle for projection out constraint components
 181      * of derivatives of the coordinates.
 182      * Berk Hess 2008-1-10
 183      */
 184
 185     settleparam_t *p;
 186     real   imO,imH,dOH,dHH,invdOH,invdHH;
 187     matrix invmat;
 188     int    i,m,m2,ow1,hw2,hw3;
 189     rvec   roh2,roh3,rhh,dc,fc,fcv;
 190     rvec   derm[3],derpm[3];
 191     real   invvscale,vscale_nhc,veta;
 192     real   kfacOH,kfacHH;
 193
 194     calcvir_atom_end *= DIM;
 195
 196     if (econq == econqForce)
 197     {
 198         p = &settled->mass1;
 199     }
 200     else
 201     {
 202         p = &settled->massw;
 203     }
 204     imO    = p->imO;
 205     imH    = p->imH;
 206     copy_mat(p->invmat,invmat);
 207     dOH    = p->dOH;
 208     dHH    = p->dHH;
 209     invdOH = p->invdOH;
 210     invdHH = p->invdHH;
 211
 212     veta = vetavar->veta;
 213     vscale_nhc = vetavar->vscale_nhc[0]; /* assume the first temperature control group. */
 214
 215 #ifdef PRAGMAS
 216 #pragma ivdep
 217 #endif
 218
 219     for (i=0; i<nsettle; i++)
 220     {
 221         ow1 = iatoms[i*4+1];
 222         hw2 = iatoms[i*4+2];
 223         hw3 = iatoms[i*4+3];
 224
 225
 226         for(m=0; m<DIM; m++)
 227         {
 228             /* in the velocity case, these are the velocities, so we
 229                need to modify with the pressure control velocities! */
 230
 231             derm[0][m]  = vscale_nhc*der[ow1][m] + veta*x[ow1][m];
 232             derm[1][m]  = vscale_nhc*der[hw2][m] + veta*x[hw2][m];
 233             derm[2][m]  = vscale_nhc*der[hw3][m] + veta*x[hw3][m];
 234
 235         }
 236         /* 27 flops */
 237
 238         if (pbc == NULL)
 239         {
 240             rvec_sub(x[ow1],x[hw2],roh2);
 241             rvec_sub(x[ow1],x[hw3],roh3);
 242             rvec_sub(x[hw2],x[hw3],rhh);
 243         }
 244         else
 245         {
 246             pbc_dx_aiuc(pbc,x[ow1],x[hw2],roh2);
 247             pbc_dx_aiuc(pbc,x[ow1],x[hw3],roh3);
 248             pbc_dx_aiuc(pbc,x[hw2],x[hw3],rhh);
 249         }
 250         svmul(invdOH,roh2,roh2);
 251         svmul(invdOH,roh3,roh3);
 252         svmul(invdHH,rhh,rhh);
 253         /* 18 flops */
 254
 255         /* Determine the projections of der(modified) on the bonds */
 256         clear_rvec(dc);
 257         for(m=0; m<DIM; m++)
 258         {
 259             dc[0] += (derm[0][m] - derm[1][m])*roh2[m];
 260             dc[1] += (derm[0][m] - derm[2][m])*roh3[m];
 261             dc[2] += (derm[1][m] - derm[2][m])*rhh [m];
 262         }
 263         /* 27 flops */
 264
 265         /* Determine the correction for the three bonds */
 266         mvmul(invmat,dc,fc);
 267         /* 15 flops */
 268
 269         /* divide velocity by vscale_nhc for determining constrained velocities, since they
 270            have not yet been multiplied */
 271         svmul(1.0/vscale_nhc,fc,fcv);
 272         /* 7? flops */
 273
 274         /* Subtract the corrections from derp */
 275         for(m=0; m<DIM; m++)
 276         {
 277             derp[ow1][m] -= imO*( fcv[0]*roh2[m] + fcv[1]*roh3[m]);
 278             derp[hw2][m] -= imH*(-fcv[0]*roh2[m] + fcv[2]*rhh [m]);
 279             derp[hw3][m] -= imH*(-fcv[1]*roh3[m] - fcv[2]*rhh [m]);
 280         }
 281
 282         /* 45 flops */
 283
 284         if (ow1 < calcvir_atom_end)
 285         {
 286             /* Determining r \dot m der is easy,
 287              * since fc contains the mass weighted corrections for der.
 288              */
 289
 290             for(m=0; m<DIM; m++)
 291             {
 292                 for(m2=0; m2<DIM; m2++)
 293                 {
 294                     rmdder[m][m2] +=
 295                         dOH*roh2[m]*roh2[m2]*fcv[0] +
 296                         dOH*roh3[m]*roh3[m2]*fcv[1] +
 297                         dHH*rhh [m]*rhh [m2]*fcv[2];
 298                 }
 299             }
 300         }
 301     }
 302
 303     if (calcvir_atom_end > 0)
 304     {
 305         /* Correct rmdder, which will be used to calcualate the virial;
 306          * we need to use the unscaled multipliers in the virial.
 307          */
 308         msmul(rmdder,1.0/vetavar->vscale,rmdder);
 309     }
 310 }
 311
 312
 313 void csettle(gmx_settledata_t settled,
 314              int nsettle, t_iatom iatoms[],
 315              const t_pbc *pbc,
 316              real b4[], real after[],
 317              real invdt,real *v,int CalcVirAtomEnd,
 318              tensor rmdr,int *error,t_vetavars *vetavar)
 319 {
 320     /* ***************************************************************** */
 321     /*                                                               ** */
 322     /*    Subroutine : setlep - reset positions of TIP3P waters      ** */
 323     /*    Author : Shuichi Miyamoto                                  ** */
 324     /*    Date of last update : Oct. 1, 1992                         ** */
 325     /*                                                               ** */
 326     /*    Reference for the SETTLE algorithm                         ** */
 327     /*           S. Miyamoto et al., J. Comp. Chem., 13, 952 (1992). ** */
 328     /*                                                               ** */
 329     /* ***************************************************************** */
 330
 331     /* Initialized data */
 332     settleparam_t *p;
 333     real   mO,mH,mOs,mHs,invdts;
 334     /* These three weights need have double precision. Using single precision
 335      * can result in huge velocity and pressure deviations. */
 336     double wo,wh,wohh;
 337     real   ra,rb,rc,irc2,dOH,dHH;
 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     int  is;
 360
 361     *error = -1;
 362
 363     CalcVirAtomEnd *= 3;
 364
 365     p = &settled->massw;
 366     mO   = p->mO;
 367     mH   = p->mH;
 368     wo   = p->wo;
 369     wh   = p->wh;
 370     wohh = p->wohh;
 371     rc   = p->rc;
 372     ra   = p->ra;
 373     rb   = p->rb;
 374     irc2 = p->irc2;
 375     dOH  = p->dOH;
 376     dHH  = p->dHH;
 377
 378     mOs  = mO / vetavar->rvscale;
 379     mHs  = mH / vetavar->rvscale;
 380     invdts = invdt / vetavar->rscale;
 381
 382 #ifdef PRAGMAS
 383 #pragma ivdep
 384 #endif
 385   for (i = 0; i < nsettle; ++i) {
 386     bOK = TRUE;
 387     /*    --- Step1  A1' ---      */
 388     ow1 = iatoms[i*4+1] * 3;
 389     hw2 = iatoms[i*4+2] * 3;
 390     hw3 = iatoms[i*4+3] * 3;
 391     if (pbc == NULL)
 392     {
 393         xb0 = b4[hw2    ] - b4[ow1];
 394         yb0 = b4[hw2 + 1] - b4[ow1 + 1];
 395         zb0 = b4[hw2 + 2] - b4[ow1 + 2];
 396         xc0 = b4[hw3    ] - b4[ow1];
 397         yc0 = b4[hw3 + 1] - b4[ow1 + 1];
 398         zc0 = b4[hw3 + 2] - b4[ow1 + 2];
 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,dx);
 414         if (is == CENTRAL)
 415         {
 416             clear_rvec(sh_hw2);
 417         }
 418         else
 419         {
 420             sh_hw2[XX] = after[hw2    ] - (after[ow1    ] + dx[XX]);
 421             sh_hw2[YY] = after[hw2 + 1] - (after[ow1 + 1] + dx[YY]);
 422             sh_hw2[ZZ] = after[hw2 + 2] - (after[ow1 + 2] + dx[ZZ]);
 423             rvec_dec(after+hw2,sh_hw2);
 424         }
 425         is = pbc_dx_aiuc(pbc,after+hw3,after+ow1,dx);
 426         if (is == CENTRAL)
 427         {
 428             clear_rvec(sh_hw3);
 429         }
 430         else
 431         {
 432             sh_hw3[XX] = after[hw3    ] - (after[ow1    ] + dx[XX]);
 433             sh_hw3[YY] = after[hw3 + 1] - (after[ow1 + 1] + dx[YY]);
 434             sh_hw3[ZZ] = after[hw3 + 2] - (after[ow1 + 2] + dx[ZZ]);
 435             rvec_dec(after+hw3,sh_hw3);
 436         }
 437     }
 438
 439     xcom = (after[ow1    ] * wo + (after[hw2    ] + after[hw3    ]) * wh);
 440     ycom = (after[ow1 + 1] * wo + (after[hw2 + 1] + after[hw3 + 1]) * wh);
 441     zcom = (after[ow1 + 2] * wo + (after[hw2 + 2] + after[hw3 + 2]) * wh);
 442     /* 12 flops */
 443
 444     xa1 = after[ow1    ] - xcom;
 445     ya1 = after[ow1 + 1] - ycom;
 446     za1 = after[ow1 + 2] - zcom;
 447     xb1 = after[hw2    ] - xcom;
 448     yb1 = after[hw2 + 1] - ycom;
 449     zb1 = after[hw2 + 2] - zcom;
 450     xc1 = after[hw3    ] - xcom;
 451     yc1 = after[hw3 + 1] - ycom;
 452     zc1 = after[hw3 + 2] - zcom;
 453     /* 9 flops */
 454
 455     xakszd = yb0 * zc0 - zb0 * yc0;
 456     yakszd = zb0 * xc0 - xb0 * zc0;
 457     zakszd = xb0 * yc0 - yb0 * xc0;
 458     xaksxd = ya1 * zakszd - za1 * yakszd;
 459     yaksxd = za1 * xakszd - xa1 * zakszd;
 460     zaksxd = xa1 * yakszd - ya1 * xakszd;
 461     xaksyd = yakszd * zaksxd - zakszd * yaksxd;
 462     yaksyd = zakszd * xaksxd - xakszd * zaksxd;
 463     zaksyd = xakszd * yaksxd - yakszd * xaksxd;
 464     /* 27 flops */
 465
 466     axlng = gmx_invsqrt(xaksxd * xaksxd + yaksxd * yaksxd + zaksxd * zaksxd);
 467     aylng = gmx_invsqrt(xaksyd * xaksyd + yaksyd * yaksyd + zaksyd * zaksyd);
 468     azlng = gmx_invsqrt(xakszd * xakszd + yakszd * yakszd + zakszd * zakszd);
 469
 470     trns11 = xaksxd * axlng;
 471     trns21 = yaksxd * axlng;
 472     trns31 = zaksxd * axlng;
 473     trns12 = xaksyd * aylng;
 474     trns22 = yaksyd * aylng;
 475     trns32 = zaksyd * aylng;
 476     trns13 = xakszd * azlng;
 477     trns23 = yakszd * azlng;
 478     trns33 = zakszd * azlng;
 479     /* 24 flops */
 480
 481     xb0d = trns11 * xb0 + trns21 * yb0 + trns31 * zb0;
 482     yb0d = trns12 * xb0 + trns22 * yb0 + trns32 * zb0;
 483     xc0d = trns11 * xc0 + trns21 * yc0 + trns31 * zc0;
 484     yc0d = trns12 * xc0 + trns22 * yc0 + trns32 * zc0;
 485     /*
 486     xa1d = trns11 * xa1 + trns21 * ya1 + trns31 * za1;
 487     ya1d = trns12 * xa1 + trns22 * ya1 + trns32 * za1;
 488     */
 489     za1d = trns13 * xa1 + trns23 * ya1 + trns33 * za1;
 490     xb1d = trns11 * xb1 + trns21 * yb1 + trns31 * zb1;
 491     yb1d = trns12 * xb1 + trns22 * yb1 + trns32 * zb1;
 492     zb1d = trns13 * xb1 + trns23 * yb1 + trns33 * zb1;
 493     xc1d = trns11 * xc1 + trns21 * yc1 + trns31 * zc1;
 494     yc1d = trns12 * xc1 + trns22 * yc1 + trns32 * zc1;
 495     zc1d = trns13 * xc1 + trns23 * yc1 + trns33 * zc1;
 496     /* 65 flops */
 497
 498     sinphi = za1d * gmx_invsqrt(ra*ra);
 499     tmp    = 1.0 - sinphi * sinphi;
 500     if (tmp <= 0) {
 501       bOK = FALSE;
 502     } else {
 503       tmp2   = gmx_invsqrt(tmp);
 504       cosphi = tmp*tmp2;
 505       sinpsi = (zb1d - zc1d) * irc2 * tmp2;
 506       tmp2   = 1.0 - sinpsi * sinpsi;
 507       if (tmp2 <= 0) {
 508         bOK = FALSE;
 509       } else {
 510         cospsi = tmp2*gmx_invsqrt(tmp2);
 511       }
 512     }
 513     /* 46 flops */
 514
 515     if (bOK) {
 516       ya2d =  ra * cosphi;
 517       xb2d = -rc * cospsi;
 518       t1   = -rb * cosphi;
 519       t2   =  rc * sinpsi * sinphi;
 520       yb2d =  t1 - t2;
 521       yc2d =  t1 + t2;
 522       /* 7 flops */
 523
 524       /*     --- Step3  al,be,ga                      --- */
 525       alpa   = xb2d * (xb0d - xc0d) + yb0d * yb2d + yc0d * yc2d;
 526       beta   = xb2d * (yc0d - yb0d) + xb0d * yb2d + xc0d * yc2d;
 527       gama   = xb0d * yb1d - xb1d * yb0d + xc0d * yc1d - xc1d * yc0d;
 528       al2be2 = alpa * alpa + beta * beta;
 529       tmp2   = (al2be2 - gama * gama);
 530       sinthe = (alpa * gama - beta * tmp2*gmx_invsqrt(tmp2)) * gmx_invsqrt(al2be2*al2be2);
 531       /* 47 flops */
 532
 533       /*  --- Step4  A3' --- */
 534       tmp2  = 1.0 - sinthe * sinthe;
 535       costhe = tmp2*gmx_invsqrt(tmp2);
 536       xa3d = -ya2d * sinthe;
 537       ya3d = ya2d * costhe;
 538       za3d = za1d;
 539       xb3d = xb2d * costhe - yb2d * sinthe;
 540       yb3d = xb2d * sinthe + yb2d * costhe;
 541       zb3d = zb1d;
 542       xc3d = -xb2d * costhe - yc2d * sinthe;
 543       yc3d = -xb2d * sinthe + yc2d * costhe;
 544       zc3d = zc1d;
 545       /* 26 flops */
 546
 547       /*    --- Step5  A3 --- */
 548       xa3 = trns11 * xa3d + trns12 * ya3d + trns13 * za3d;
 549       ya3 = trns21 * xa3d + trns22 * ya3d + trns23 * za3d;
 550       za3 = trns31 * xa3d + trns32 * ya3d + trns33 * za3d;
 551       xb3 = trns11 * xb3d + trns12 * yb3d + trns13 * zb3d;
 552       yb3 = trns21 * xb3d + trns22 * yb3d + trns23 * zb3d;
 553       zb3 = trns31 * xb3d + trns32 * yb3d + trns33 * zb3d;
 554       xc3 = trns11 * xc3d + trns12 * yc3d + trns13 * zc3d;
 555       yc3 = trns21 * xc3d + trns22 * yc3d + trns23 * zc3d;
 556       zc3 = trns31 * xc3d + trns32 * yc3d + trns33 * zc3d;
 557       /* 45 flops */
 558       after[ow1] = xcom + xa3;
 559       after[ow1 + 1] = ycom + ya3;
 560       after[ow1 + 2] = zcom + za3;
 561       after[hw2] = xcom + xb3;
 562       after[hw2 + 1] = ycom + yb3;
 563       after[hw2 + 2] = zcom + zb3;
 564       after[hw3] = xcom + xc3;
 565       after[hw3 + 1] = ycom + yc3;
 566       after[hw3 + 2] = zcom + zc3;
 567       /* 9 flops */
 568
 569       if (pbc != NULL)
 570       {
 571           rvec_inc(after+hw2,sh_hw2);
 572           rvec_inc(after+hw3,sh_hw3);
 573       }
 574
 575       dax = xa3 - xa1;
 576       day = ya3 - ya1;
 577       daz = za3 - za1;
 578       dbx = xb3 - xb1;
 579       dby = yb3 - yb1;
 580       dbz = zb3 - zb1;
 581       dcx = xc3 - xc1;
 582       dcy = yc3 - yc1;
 583       dcz = zc3 - zc1;
 584       /* 9 flops, counted with the virial */
 585
 586       if (v != NULL) {
 587           v[ow1]     += dax*invdts;
 588           v[ow1 + 1] += day*invdts;
 589           v[ow1 + 2] += daz*invdts;
 590           v[hw2]     += dbx*invdts;
 591           v[hw2 + 1] += dby*invdts;
 592           v[hw2 + 2] += dbz*invdts;
 593           v[hw3]     += dcx*invdts;
 594           v[hw3 + 1] += dcy*invdts;
 595           v[hw3 + 2] += dcz*invdts;
 596           /* 3*6 flops */
 597       }
 598
 599       if (ow1 < CalcVirAtomEnd) {
 600           mdax = mOs*dax;
 601           mday = mOs*day;
 602           mdaz = mOs*daz;
 603           mdbx = mHs*dbx;
 604           mdby = mHs*dby;
 605           mdbz = mHs*dbz;
 606           mdcx = mHs*dcx;
 607           mdcy = mHs*dcy;
 608           mdcz = mHs*dcz;
 609           rmdr[XX][XX] -= b4[ow1  ]*mdax + (b4[ow1  ]+xb0)*mdbx + (b4[ow1  ]+xc0)*mdcx;
 610           rmdr[XX][YY] -= b4[ow1  ]*mday + (b4[ow1  ]+xb0)*mdby + (b4[ow1  ]+xc0)*mdcy;
 611           rmdr[XX][ZZ] -= b4[ow1  ]*mdaz + (b4[ow1  ]+xb0)*mdbz + (b4[ow1  ]+xc0)*mdcz;
 612           rmdr[YY][XX] -= b4[ow1+1]*mdax + (b4[ow1+1]+yb0)*mdbx + (b4[ow1+1]+yc0)*mdcx;
 613           rmdr[YY][YY] -= b4[ow1+1]*mday + (b4[ow1+1]+yb0)*mdby + (b4[ow1+1]+yc0)*mdcy;
 614           rmdr[YY][ZZ] -= b4[ow1+1]*mdaz + (b4[ow1+1]+yb0)*mdbz + (b4[ow1+1]+yc0)*mdcz;
 615           rmdr[ZZ][XX] -= b4[ow1+2]*mdax + (b4[ow1+2]+zb0)*mdbx + (b4[ow1+2]+zc0)*mdcx;
 616           rmdr[ZZ][YY] -= b4[ow1+2]*mday + (b4[ow1+2]+zb0)*mdby + (b4[ow1+2]+zc0)*mdcy;
 617           rmdr[ZZ][ZZ] -= b4[ow1+2]*mdaz + (b4[ow1+2]+zb0)*mdbz + (b4[ow1+2]+zc0)*mdcz;
 618         /* 3*24 - 9 flops */
 619       }
 620     } else {
 621       *error = i;
 622     }
 623 #ifdef DEBUG
 624     if (debug)
 625     {
 626         check_cons(debug,"settle",after,ow1,hw2,hw3);
 627     }
 628 #endif
 629   }
 630 }