include/vec.h

   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  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29
  30 #ifndef _vec_h
  31 #define _vec_h
  32
  33 static char *SRCID_vec_h = "$Id$";
  34
  35 #ifdef HAVE_IDENT
  36 #ident  "@(#) vec.h 1.8 12/16/92"
  37 #endif /* HAVE_IDENT */
  38
  39 #include "maths.h"
  40 #include "typedefs.h"
  41 #include "sysstuff.h"
  42 #include "macros.h"
  43 #include "fatal.h"
  44 #include "lutab.h"
  45
  46 #ifdef _lnx_
  47 #define gmx_inline inline
  48 #else
  49 #define gmx_inline
  50 #endif
  51
  52 #ifdef CINVSQRT
  53 static gmx_inline real invsqrt(float x)
  54 {
  55   const real  half=0.5;
  56   const real  three=3.0;
  57   t_convert   result,bit_pattern;
  58   word        exp,fract;
  59   float       lu;
  60   real        y;
  61 #ifdef DOUBLE
  62   real        y2;
  63 #endif
  64
  65   bit_pattern.fval=x;
  66   exp   = EXP_ADDR(bit_pattern.bval);
  67   fract = FRACT_ADDR(bit_pattern.bval);
  68   result.bval=lookup_table.exp_seed[exp] | lookup_table.fract_seed[fract];
  69   lu    = result.fval;
  70
  71   y=(half*lu*(three-((x*lu)*lu)));
  72 #ifdef DOUBLE
  73   y2=(half*y*(three-((x*y)*y)));
  74
  75   return y2;                    /* 10 Flops */
  76 #else
  77   return y;                     /* 5  Flops */
  78 #endif
  79 }
  80 #define INVSQRT_DONE
  81 #endif
  82
  83 #ifndef INVSQRT_DONE
  84 #define invsqrt(x) (1.0f/sqrt(x))
  85 #endif
  86
  87 static gmx_inline real sqr(real x)
  88 {
  89   return (x*x);
  90 }
  91
  92 static gmx_inline void rvec_add(rvec a,rvec b,rvec c)
  93 {
  94   real x,y,z;
  95
  96   x=a[XX]+b[XX];
  97   y=a[YY]+b[YY];
  98   z=a[ZZ]+b[ZZ];
  99
 100   c[XX]=x;
 101   c[YY]=y;
 102   c[ZZ]=z;
 103 }
 104
 105 static gmx_inline void rvec_inc(rvec a,rvec b)
 106 {
 107   real x,y,z;
 108
 109   x=a[XX]+b[XX];
 110   y=a[YY]+b[YY];
 111   z=a[ZZ]+b[ZZ];
 112
 113   a[XX]=x;
 114   a[YY]=y;
 115   a[ZZ]=z;
 116 }
 117
 118 static gmx_inline void rvec_sub(rvec a,rvec b,rvec c)
 119 {
 120   real x,y,z;
 121
 122   x=a[XX]-b[XX];
 123   y=a[YY]-b[YY];
 124   z=a[ZZ]-b[ZZ];
 125
 126   c[XX]=x;
 127   c[YY]=y;
 128   c[ZZ]=z;
 129 }
 130
 131 static gmx_inline void rvec_dec(rvec a,rvec b)
 132 {
 133   real x,y,z;
 134
 135   x=a[XX]-b[XX];
 136   y=a[YY]-b[YY];
 137   z=a[ZZ]-b[ZZ];
 138
 139   a[XX]=x;
 140   a[YY]=y;
 141   a[ZZ]=z;
 142 }
 143
 144 static gmx_inline void copy_rvec(rvec a,rvec b)
 145 {
 146   b[XX]=a[XX];
 147   b[YY]=a[YY];
 148   b[ZZ]=a[ZZ];
 149 }
 150
 151 static gmx_inline void copy_mat(matrix a,matrix b)
 152 {
 153   copy_rvec(a[XX],b[XX]);
 154   copy_rvec(a[YY],b[YY]);
 155   copy_rvec(a[ZZ],b[ZZ]);
 156 }
 157
 158 static gmx_inline void svmul(real a,rvec v1,rvec v2)
 159 {
 160   v2[XX]=a*v1[XX];
 161   v2[YY]=a*v1[YY];
 162   v2[ZZ]=a*v1[ZZ];
 163 }
 164
 165 static gmx_inline real distance2(rvec v1, rvec v2)
 166 {
 167   return sqr(v2[XX]-v1[XX]) + sqr(v2[YY]-v1[YY]) + sqr(v2[ZZ]-v1[ZZ]);
 168 }
 169
 170 static gmx_inline void clear_rvec(rvec a)
 171 {
 172   const real nul=0.0;
 173
 174   a[XX]=a[YY]=a[ZZ]=nul;
 175 }
 176
 177 static gmx_inline void clear_rvecs(int n,rvec v[])
 178 {
 179   int i;
 180
 181   for(i=0; (i<n); i++)
 182     clear_rvec(v[i]);
 183 }
 184
 185 static gmx_inline void clear_mat(matrix a)
 186 {
 187   const real nul=0.0;
 188
 189   a[XX][XX]=a[XX][YY]=a[XX][ZZ]=nul;
 190   a[YY][XX]=a[YY][YY]=a[YY][ZZ]=nul;
 191   a[ZZ][XX]=a[ZZ][YY]=a[ZZ][ZZ]=nul;
 192 }
 193
 194 static gmx_inline real iprod(rvec a,rvec b)
 195 {
 196   return (a[XX]*b[XX]+a[YY]*b[YY]+a[ZZ]*b[ZZ]);
 197 }
 198
 199 static gmx_inline real norm2(rvec a)
 200 {
 201   return a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ];
 202 }
 203
 204 static gmx_inline real norm(rvec a)
 205 {
 206   return sqrt(a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ]);
 207 }
 208
 209 static gmx_inline real cos_angle(rvec a,rvec b)
 210 {
 211   /*
 212    *                  ax*bx + ay*by + az*bz
 213    * cos-vec (a,b) =  ---------------------
 214    *                      ||a|| * ||b||
 215    */
 216   real   cos;
 217   int    m;
 218   double aa,bb,ip,ipa,ipb; /* For accuracy these must be double! */
 219
 220   ip=ipa=ipb=0;
 221   for(m=0; (m<DIM); m++) {              /* 18           */
 222     aa   = a[m];
 223     bb   = b[m];
 224     ip  += aa*bb;
 225     ipa += aa*aa;
 226     ipb += bb*bb;
 227   }
 228   cos=ip*invsqrt(ipa*ipb);              /*  7           */
 229                                         /* 25 TOTAL     */
 230   if (cos > 1.0)
 231     return  1.0;
 232   if (cos <-1.0)
 233     return -1.0;
 234
 235   return cos;
 236 }
 237
 238 static gmx_inline void oprod(rvec a,rvec b,rvec c)
 239 {
 240   c[XX]=a[YY]*b[ZZ]-a[ZZ]*b[YY];
 241   c[YY]=a[ZZ]*b[XX]-a[XX]*b[ZZ];
 242   c[ZZ]=a[XX]*b[YY]-a[YY]*b[XX];
 243 }
 244
 245 static gmx_inline void mmul(matrix a,matrix b,matrix dest)
 246 {
 247   dest[XX][XX]=a[XX][XX]*b[XX][XX]+a[XX][YY]*b[YY][XX]+a[XX][ZZ]*b[ZZ][XX];
 248   dest[YY][XX]=a[YY][XX]*b[XX][XX]+a[YY][YY]*b[YY][XX]+a[YY][ZZ]*b[ZZ][XX];
 249   dest[ZZ][XX]=a[ZZ][XX]*b[XX][XX]+a[ZZ][YY]*b[YY][XX]+a[ZZ][ZZ]*b[ZZ][XX];
 250   dest[XX][YY]=a[XX][XX]*b[XX][YY]+a[XX][YY]*b[YY][YY]+a[XX][ZZ]*b[ZZ][YY];
 251   dest[YY][YY]=a[YY][XX]*b[XX][YY]+a[YY][YY]*b[YY][YY]+a[YY][ZZ]*b[ZZ][YY];
 252   dest[ZZ][YY]=a[ZZ][XX]*b[XX][YY]+a[ZZ][YY]*b[YY][YY]+a[ZZ][ZZ]*b[ZZ][YY];
 253   dest[XX][ZZ]=a[XX][XX]*b[XX][ZZ]+a[XX][YY]*b[YY][ZZ]+a[XX][ZZ]*b[ZZ][ZZ];
 254   dest[YY][ZZ]=a[YY][XX]*b[XX][ZZ]+a[YY][YY]*b[YY][ZZ]+a[YY][ZZ]*b[ZZ][ZZ];
 255   dest[ZZ][ZZ]=a[ZZ][XX]*b[XX][ZZ]+a[ZZ][YY]*b[YY][ZZ]+a[ZZ][ZZ]*b[ZZ][ZZ];
 256 }
 257
 258 static gmx_inline real det(matrix a)
 259 {
 260   return ( a[XX][XX]*(a[YY][YY]*a[ZZ][ZZ]-a[ZZ][YY]*a[YY][ZZ])
 261           -a[YY][XX]*(a[XX][YY]*a[ZZ][ZZ]-a[ZZ][YY]*a[XX][ZZ])
 262           +a[ZZ][XX]*(a[XX][YY]*a[YY][ZZ]-a[YY][YY]*a[XX][ZZ]));
 263 }
 264
 265 static gmx_inline void m_add(matrix a,matrix b,matrix dest)
 266 {
 267   dest[XX][XX]=a[XX][XX]+b[XX][XX];
 268   dest[XX][YY]=a[XX][YY]+b[XX][YY];
 269   dest[XX][ZZ]=a[XX][ZZ]+b[XX][ZZ];
 270   dest[YY][XX]=a[YY][XX]+b[YY][XX];
 271   dest[YY][YY]=a[YY][YY]+b[YY][YY];
 272   dest[YY][ZZ]=a[YY][ZZ]+b[YY][ZZ];
 273   dest[ZZ][XX]=a[ZZ][XX]+b[ZZ][XX];
 274   dest[ZZ][YY]=a[ZZ][YY]+b[ZZ][YY];
 275   dest[ZZ][ZZ]=a[ZZ][ZZ]+b[ZZ][ZZ];
 276 }
 277
 278 static gmx_inline void m_sub(matrix a,matrix b,matrix dest)
 279 {
 280   dest[XX][XX]=a[XX][XX]-b[XX][XX];
 281   dest[XX][YY]=a[XX][YY]-b[XX][YY];
 282   dest[XX][ZZ]=a[XX][ZZ]-b[XX][ZZ];
 283   dest[YY][XX]=a[YY][XX]-b[YY][XX];
 284   dest[YY][YY]=a[YY][YY]-b[YY][YY];
 285   dest[YY][ZZ]=a[YY][ZZ]-b[YY][ZZ];
 286   dest[ZZ][XX]=a[ZZ][XX]-b[ZZ][XX];
 287   dest[ZZ][YY]=a[ZZ][YY]-b[ZZ][YY];
 288   dest[ZZ][ZZ]=a[ZZ][ZZ]-b[ZZ][ZZ];
 289 }
 290
 291 static gmx_inline void msmul(matrix m1,real r1,matrix dest)
 292 {
 293   dest[XX][XX]=r1*m1[XX][XX];
 294   dest[XX][YY]=r1*m1[XX][YY];
 295   dest[XX][ZZ]=r1*m1[XX][ZZ];
 296   dest[YY][XX]=r1*m1[YY][XX];
 297   dest[YY][YY]=r1*m1[YY][YY];
 298   dest[YY][ZZ]=r1*m1[YY][ZZ];
 299   dest[ZZ][XX]=r1*m1[ZZ][XX];
 300   dest[ZZ][YY]=r1*m1[ZZ][YY];
 301   dest[ZZ][ZZ]=r1*m1[ZZ][ZZ];
 302 }
 303
 304 static gmx_inline void m_inv(matrix src,matrix dest)
 305 {
 306   const real smallreal = 1.0e-18;
 307   const real largereal = 1.0e18;
 308   real  deter,c,fc;
 309
 310   deter=det(src);
 311   c=1.0/deter;
 312   fc=fabs(c);
 313
 314   if ((fc <= smallreal) || (fc >= largereal))
 315     fatal_error(0,"Determinant = %f",1.0/c);
 316
 317   dest[XX][XX]= c*(src[YY][YY]*src[ZZ][ZZ]-src[ZZ][YY]*src[YY][ZZ]);
 318   dest[XX][YY]=-c*(src[XX][YY]*src[ZZ][ZZ]-src[ZZ][YY]*src[XX][ZZ]);
 319   dest[XX][ZZ]= c*(src[XX][YY]*src[YY][ZZ]-src[YY][YY]*src[XX][ZZ]);
 320   dest[YY][XX]=-c*(src[YY][XX]*src[ZZ][ZZ]-src[ZZ][XX]*src[YY][ZZ]);
 321   dest[YY][YY]= c*(src[XX][XX]*src[ZZ][ZZ]-src[ZZ][XX]*src[XX][ZZ]);
 322   dest[YY][ZZ]=-c*(src[XX][XX]*src[YY][ZZ]-src[YY][XX]*src[XX][ZZ]);
 323   dest[ZZ][XX]= c*(src[YY][XX]*src[ZZ][YY]-src[ZZ][XX]*src[YY][YY]);
 324   dest[ZZ][YY]=-c*(src[XX][XX]*src[ZZ][YY]-src[ZZ][XX]*src[XX][YY]);
 325   dest[ZZ][ZZ]= c*(src[XX][XX]*src[YY][YY]-src[YY][XX]*src[XX][YY]);
 326 }
 327
 328 static gmx_inline void mvmul(matrix a,rvec src,rvec dest)
 329 {
 330   dest[XX]=a[XX][XX]*src[XX]+a[XX][YY]*src[YY]+a[XX][ZZ]*src[ZZ];
 331   dest[YY]=a[YY][XX]*src[XX]+a[YY][YY]*src[YY]+a[YY][ZZ]*src[ZZ];
 332   dest[ZZ]=a[ZZ][XX]*src[XX]+a[ZZ][YY]*src[YY]+a[ZZ][ZZ]*src[ZZ];
 333 }
 334
 335 static gmx_inline void unitv(rvec src,rvec dest)
 336 {
 337   real linv;
 338
 339   linv=invsqrt(iprod(src,src));
 340   dest[XX]=linv*src[XX];
 341   dest[YY]=linv*src[YY];
 342   dest[ZZ]=linv*src[ZZ];
 343 }
 344
 345 static gmx_inline real trace(matrix m)
 346 {
 347   return (m[XX][XX]+m[YY][YY]+m[ZZ][ZZ]);
 348 }
 349
 350 static gmx_inline real _divide(real a,real b,char *file,int line)
 351 {
 352   if (b == 0.0)
 353     fatal_error(0,"Dividing by zero, file %s, line %d",file,line);
 354   return a/b;
 355 }
 356
 357 static gmx_inline int _mod(int a,int b,char *file,int line)
 358 {
 359   if (b == 0)
 360     fatal_error(0,"Modulo zero, file %s, line %d",file,line);
 361   return a % b;
 362 }
 363
 364 #define divide(a,b) _divide((a),(b),__FILE__,__LINE__)
 365 #define mod(a,b)    _mod((a),(b),__FILE__,__LINE__)
 366
 367 #endif  /* _vec_h */
 368
 369