src/mdlib/cshake.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_cshake_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "sysstuff.h"
  33 #include "typedefs.h"
  34 #include "smalloc.h"
  35 #include "pbc.h"
  36 #include "txtdump.h"
  37 #include "constr.h"
  38 #include "vec.h"
  39 #include "nrnb.h"
  40
  41 void cshake(atom_id iatom[],int ncon,int *nnit,int maxnit,
  42             real dist2[],real xp[],real rij[],real m2[],
  43             real invmass[],real tt[],int *nerror)
  44 {
  45   const   real mytol=1e-6;
  46
  47   int     ll,i,j,i3,j3,l3;
  48   int     ix,iy,iz,jx,jy,jz;
  49   real    toler,rpij2,rrpr,tx,ty,tz,diff,acor,im,jm;
  50   real    xh,yh,zh,rijx,rijy,rijz;
  51   real    tix,tiy,tiz;
  52   real    tjx,tjy,tjz;
  53   int     nit,error,iconv,nconv;
  54
  55   error=0;
  56   nconv=1;
  57   for (nit=0; (nit<maxnit) && (nconv != 0) && (error == 0); nit++) {
  58     nconv=0;
  59     for(ll=0; (ll<ncon) && (error == 0); ll++) {
  60       l3    = 3*ll;
  61       rijx  = rij[l3+XX];
  62       rijy  = rij[l3+YY];
  63       rijz  = rij[l3+ZZ];
  64       i     = iatom[l3+1];
  65       j     = iatom[l3+2];
  66       i3    = 3*i;
  67       j3    = 3*j;
  68       ix    = i3+XX;
  69       iy    = i3+YY;
  70       iz    = i3+ZZ;
  71       jx    = j3+XX;
  72       jy    = j3+YY;
  73       jz    = j3+ZZ;
  74
  75       tx      = xp[ix]-xp[jx];
  76       ty      = xp[iy]-xp[jy];
  77       tz      = xp[iz]-xp[jz];
  78       rpij2   = tx*tx+ty*ty+tz*tz;
  79       toler   = dist2[ll];
  80       diff    = toler-rpij2;
  81
  82       /* iconv is zero when the error is smaller than a bound */
  83       iconv   = fabs(diff)*tt[ll];
  84
  85       if (iconv != 0) {
  86         nconv   = nconv + iconv;
  87         rrpr    = rijx*tx+rijy*ty+rijz*tz;
  88
  89         if (rrpr < toler*mytol)
  90           error=ll;
  91         else {
  92           acor    = diff*m2[ll]/rrpr;
  93           im      = invmass[i];
  94           jm      = invmass[j];
  95           xh      = rijx*acor;
  96           yh      = rijy*acor;
  97           zh      = rijz*acor;
  98           tix     = xp[ix] + xh*im;
  99           tiy     = xp[iy] + yh*im;
 100           tiz     = xp[iz] + zh*im;
 101           tjx     = xp[jx] - xh*jm;
 102           tjy     = xp[jy] - yh*jm;
 103           tjz     = xp[jz] - zh*jm;
 104           xp[ix] = tix;
 105           xp[iy] = tiy;
 106           xp[iz] = tiz;
 107           xp[jx] = tjx;
 108           xp[jy] = tjy;
 109           xp[jz] = tjz;
 110         }
 111       }
 112     }
 113   }
 114   *nnit=nit;
 115   *nerror=error;
 116 }
 117