src/mdlib/clincs.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_clincs_c = "$Id$";
  30
  31 #include <math.h>
  32 #include "main.h"
  33 #include "constr.h"
  34 #include "physics.h"
  35 #include "vec.h"
  36
  37 void clincs(rvec *x,rvec *xp,int ncons,
  38             int *bla1,int *bla2,int *blnr,int *blbnb,real *bllen,
  39             real *blc,real *blcc,real *blm,
  40             int nit,int nrec,real *invmass,rvec * r,
  41             real *rhs1,real *rhs2,real *sol,real wangle,int *warn,
  42             real *lambda)
  43 {
  44   int     b,i,j,k,n,it,rec;
  45   real    tmp0,tmp1,tmp2,im1,im2,mvb,rlen,len,wfac,lam;
  46   real    u0,u1,u2,v0,v1,v2;
  47   real    *tmp;
  48
  49   *warn=0;
  50
  51   /* Compute normalized i-j vectors */
  52   for(b=0;b<ncons;b++) {
  53     i=bla1[b];
  54     j=bla2[b];
  55     tmp0=x[i][0]-x[j][0];
  56     tmp1=x[i][1]-x[j][1];
  57     tmp2=x[i][2]-x[j][2];
  58     rlen=invsqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
  59     r[b][0]=rlen*tmp0;
  60     r[b][1]=rlen*tmp1;
  61     r[b][2]=rlen*tmp2;
  62   } /* 16 ncons flops */
  63
  64   for(b=0;b<ncons;b++) {
  65     tmp0=r[b][0];
  66     tmp1=r[b][1];
  67     tmp2=r[b][2];
  68     len=bllen[b];
  69     i=bla1[b];
  70     j=bla2[b];
  71     for(n=blnr[b];n<blnr[b+1];n++) {
  72       k=blbnb[n];
  73       blm[n]=blcc[n]*(tmp0*r[k][0]+tmp1*r[k][1]+tmp2*r[k][2]);
  74     } /* 6 nr flops */
  75     mvb=blc[b]*(tmp0*(xp[i][0]-xp[j][0])+
  76                 tmp1*(xp[i][1]-xp[j][1])+
  77                 tmp2*(xp[i][2]-xp[j][2])-len);
  78     rhs1[b]=mvb;
  79     sol[b]=mvb;
  80     /* 8 flops */
  81   }
  82   /* Together: 24*ncons + 6*nrtot flops */
  83
  84
  85   for(rec=0;rec<nrec;rec++) {
  86     for(b=0;b<ncons;b++) {
  87       mvb=0;
  88       for(n=blnr[b];n<blnr[b+1];n++) {
  89         j=blbnb[n];
  90         mvb=mvb+blm[n]*rhs1[j];
  91       }
  92       rhs2[b]=mvb;
  93       sol[b]=sol[b]+mvb;
  94     }
  95     tmp=rhs1;
  96     rhs1=rhs2;
  97     rhs2=tmp;
  98   } /* nrec*(ncons+2*nrtot) flops */
  99
 100   for(b=0;b<ncons;b++) {
 101     i=bla1[b];
 102     j=bla2[b];
 103     mvb=blc[b]*sol[b];
 104     lambda[b]=-mvb;
 105     im1=invmass[i];
 106     im2=invmass[j];
 107     tmp0=r[b][0]*mvb;
 108     tmp1=r[b][1]*mvb;
 109     tmp2=r[b][2]*mvb;
 110     u0=xp[i][0]-tmp0*im1;
 111     u1=xp[i][1]-tmp1*im1;
 112     u2=xp[i][2]-tmp2*im1;
 113     v0=xp[j][0]+tmp0*im2;
 114     v1=xp[j][1]+tmp1*im2;
 115     v2=xp[j][2]+tmp2*im2;
 116     xp[i][0]=u0;
 117     xp[i][1]=u1;
 118     xp[i][2]=u2;
 119     xp[j][0]=v0;
 120     xp[j][1]=v1;
 121     xp[j][2]=v2;
 122   } /* 16 ncons flops */
 123
 124
 125
 126   /*
 127   ********  Correction for centripetal effects  ********
 128   */
 129
 130   wfac=cos(DEG2RAD*wangle);
 131   wfac=wfac*wfac;
 132
 133   for(it=0; it<nit; it++) {
 134
 135     for(b=0;b<ncons;b++) {
 136       len=bllen[b];
 137       i=bla1[b];
 138       j=bla2[b];
 139       tmp0=xp[i][0]-xp[j][0];
 140       tmp1=xp[i][1]-xp[j][1];
 141       tmp2=xp[i][2]-xp[j][2];
 142       u1=len*len;
 143       u0=2.*u1-(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
 144       if (u0 < wfac*u1) *warn=b;
 145       if (u0 < 0) u0=0;
 146       mvb=blc[b]*(len-sqrt(u0));
 147       rhs1[b]=mvb;
 148       sol[b]=mvb;
 149     } /* 18*ncons flops */
 150
 151     for(rec=0;rec<nrec;rec++) {
 152       for(b=0;b<ncons;b++) {
 153         mvb=0;
 154         for(n=blnr[b];n<blnr[b+1];n++) {
 155           j=blbnb[n];
 156           mvb=mvb+blm[n]*rhs1[j];
 157         }
 158         rhs2[b]=mvb;
 159         sol[b]=sol[b]+mvb;
 160       }
 161       tmp=rhs1;
 162       rhs1=rhs2;
 163       rhs2=tmp;
 164     } /* nrec*(ncons+2*nrtot) flops */
 165
 166     for(b=0;b<ncons;b++) {
 167       i=bla1[b];
 168       j=bla2[b];
 169       lam=lambda[b];
 170       mvb=blc[b]*sol[b];
 171       lambda[b]=lam-mvb;
 172       im1=invmass[i];
 173       im2=invmass[j];
 174       tmp0=r[b][0]*mvb;
 175       tmp1=r[b][1]*mvb;
 176       tmp2=r[b][2]*mvb;
 177       u0=xp[i][0]-tmp0*im1;
 178       u1=xp[i][1]-tmp1*im1;
 179       u2=xp[i][2]-tmp2*im1;
 180       v0=xp[j][0]+tmp0*im2;
 181       v1=xp[j][1]+tmp1*im2;
 182       v2=xp[j][2]+tmp2*im2;
 183       xp[i][0]=u0;
 184       xp[i][1]=u1;
 185       xp[i][2]=u2;
 186       xp[j][0]=v0;
 187       xp[j][1]=v1;
 188       xp[j][2]=v2;
 189     } /* 17 ncons flops */
 190   } /* nit*ncons*(35+9*nrec) flops */
 191   /* Total:
 192    * 24*ncons + 6*nrtot + nrec*(ncons+2*nrtot)
 193    * + nit * (18*ncons + nrec*(ncons+2*nrtot) + 17 ncons)
 194    *
 195    * (24+nrec)*ncons + (6+2*nrec)*nrtot
 196    * + nit * ((35+nrec)*ncons + 2*nrec*nrtot)
 197    * if nit=1
 198    * (59+nrec)*ncons + (6+4*nrec)*nrtot
 199    */
 200 }
 201
 202 void cconerr(real *max,real *rms,int *imax,rvec *xprime,
 203              int ncons,int *bla1,int *bla2,real *bllen)
 204
 205 {
 206   real      len,d,ma,ms,tmp0,tmp1,tmp2;
 207   int       b,i,j,im;
 208
 209   ma=0;
 210   ms=0;
 211   im=0;
 212   for(b=0;b<ncons;b++) {
 213     i=bla1[b];
 214     j=bla2[b];
 215     tmp0=xprime[i][0]-xprime[j][0];
 216     tmp1=xprime[i][1]-xprime[j][1];
 217     tmp2=xprime[i][2]-xprime[j][2];
 218     len=sqrt(tmp0*tmp0+tmp1*tmp1+tmp2*tmp2);
 219     d=fabs(len/bllen[b]-1);
 220     if (d > ma) {
 221       ma=d;
 222       im=b;
 223     }
 224     ms=ms+d*d;
 225   }
 226   *max=ma;
 227   *rms=sqrt(ms/ncons);
 228   *imax=im;
 229 }
 230
 231 void lincs_warning(rvec *x,rvec *xprime,
 232                    int ncons,int *bla1,int *bla2,real *bllen,real wangle)
 233 {
 234   int b,i,j;
 235   rvec v0,v1;
 236   real wfac,d0,d1,cosine;
 237   char buf[STRLEN];
 238
 239   wfac=cos(DEG2RAD*wangle);
 240
 241   sprintf(buf,"bonds that rotated more than %g degrees:\n"
 242           " atom 1 atom 2  angle  previous, current, constraint length\n",
 243           wangle);
 244   fprintf(stderr,buf);
 245   fprintf(stdlog,buf);
 246
 247   for(b=0;b<ncons;b++) {
 248     i=bla1[b];
 249     j=bla2[b];
 250     rvec_sub(x[i],x[j],v0);
 251     rvec_sub(xprime[i],xprime[j],v1);
 252     d0=norm(v0);
 253     d1=norm(v1);
 254     cosine=iprod(v0,v1)/(d0*d1);
 255     if (cosine<wfac) {
 256       sprintf(buf," %6d %6d  %5.1f  %8.4f %8.4f    %8.4f\n",
 257               i+1,j+1,RAD2DEG*acos(cosine),d0,d1,bllen[b]);
 258       fprintf(stderr,buf);
 259       fprintf(stdlog,buf);
 260     }
 261   }
 262 }
 263
 264