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