src/gromacs/legacyheaders/vec.h

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gromacs Runs On Most of All Computer Systems
  34  */
  35 #ifndef _vec_h
  36 #define _vec_h
  37
  38 /*
  39   collection of in-line ready operations:
  40
  41   lookup-table optimized scalar operations:
  42   real gmx_invsqrt(real x)
  43   void vecinvsqrt(real in[],real out[],int n)
  44   void vecrecip(real in[],real out[],int n)
  45   real sqr(real x)
  46   double dsqr(double x)
  47
  48   vector operations:
  49   void rvec_add(const rvec a,const rvec b,rvec c)  c = a + b
  50   void dvec_add(const dvec a,const dvec b,dvec c)  c = a + b
  51   void ivec_add(const ivec a,const ivec b,ivec c)  c = a + b
  52   void rvec_inc(rvec a,const rvec b)               a += b
  53   void dvec_inc(dvec a,const dvec b)               a += b
  54   void ivec_inc(ivec a,const ivec b)               a += b
  55   void rvec_sub(const rvec a,const rvec b,rvec c)  c = a - b
  56   void dvec_sub(const dvec a,const dvec b,dvec c)  c = a - b
  57   void rvec_dec(rvec a,rvec b)                     a -= b
  58   void copy_rvec(const rvec a,rvec b)              b = a (reals)
  59   void copy_dvec(const dvec a,dvec b)              b = a (reals)
  60   void copy_ivec(const ivec a,ivec b)              b = a (integers)
  61   void ivec_sub(const ivec a,const ivec b,ivec c)  c = a - b
  62   void svmul(real a,rvec v1,rvec v2)               v2 = a * v1
  63   void dsvmul(double a,dvec v1,dvec v2)            v2 = a * v1
  64   void clear_rvec(rvec a)                          a = 0
  65   void clear_dvec(dvec a)                          a = 0
  66   void clear_ivec(rvec a)                          a = 0
  67   void clear_rvecs(int n,rvec v[])
  68   real iprod(rvec a,rvec b)                        = a . b (inner product)
  69   double diprod(dvec a,dvec b)                     = a . b (inner product)
  70   real iiprod(ivec a,ivec b)                       = a . b (integers)
  71   real norm2(rvec a)                               = | a |^2 ( = x*y*z )
  72   double dnorm2(dvec a)                            = | a |^2 ( = x*y*z )
  73   real norm(rvec a)                                = | a |
  74   double dnorm(dvec a)                             = | a |
  75   void cprod(rvec a,rvec b,rvec c)                 c = a x b (cross product)
  76   void dprod(rvec a,rvec b,rvec c)                 c = a x b (cross product)
  77   void dprod(rvec a,rvec b,rvec c)                 c = a * b (direct product)
  78   real cos_angle(rvec a,rvec b)
  79   real cos_angle_no_table(rvec a,rvec b)
  80   real distance2(rvec v1, rvec v2)                 = | v2 - v1 |^2
  81   void unitv(rvec src,rvec dest)                   dest = src / |src|
  82   void unitv_no_table(rvec src,rvec dest)          dest = src / |src|
  83
  84   matrix (3x3) operations:
  85     ! indicates that dest should not be the same as a, b or src
  86     the _ur0 varieties work on matrices that have only zeros
  87     in the upper right part, such as box matrices, these varieties
  88     could produce less rounding errors, not due to the operations themselves,
  89     but because the compiler can easier recombine the operations
  90   void copy_mat(matrix a,matrix b)                 b = a
  91   void clear_mat(matrix a)                         a = 0
  92   void mmul(matrix a,matrix b,matrix dest)      !  dest = a . b
  93   void mmul_ur0(matrix a,matrix b,matrix dest)     dest = a . b
  94   void transpose(matrix src,matrix dest)        !  dest = src*
  95   void tmmul(matrix a,matrix b,matrix dest)     !  dest = a* . b
  96   void mtmul(matrix a,matrix b,matrix dest)     !  dest = a . b*
  97   real det(matrix a)                               = det(a)
  98   void m_add(matrix a,matrix b,matrix dest)        dest = a + b
  99   void m_sub(matrix a,matrix b,matrix dest)        dest = a - b
 100   void msmul(matrix m1,real r1,matrix dest)        dest = r1 * m1
 101   void m_inv_ur0(matrix src,matrix dest)           dest = src^-1
 102   void m_inv(matrix src,matrix dest)            !  dest = src^-1
 103   void mvmul(matrix a,rvec src,rvec dest)       !  dest = a . src
 104   void mvmul_ur0(matrix a,rvec src,rvec dest)      dest = a . src
 105   void tmvmul_ur0(matrix a,rvec src,rvec dest)     dest = a* . src
 106   real trace(matrix m)                             = trace(m)
 107 */
 108
 109 #include "types/simple.h"
 110 #include "maths.h"
 111 #include "typedefs.h"
 112 #include "sysstuff.h"
 113 #include "gmx_fatal.h"
 114 #include "physics.h"
 115
 116 #ifdef __cplusplus
 117 extern "C" {
 118 #elif 0
 119 } /* avoid screwing up indentation */
 120 #endif
 121
 122
 123 #define EXP_LSB         0x00800000
 124 #define EXP_MASK        0x7f800000
 125 #define EXP_SHIFT       23
 126 #define FRACT_MASK      0x007fffff
 127 #define FRACT_SIZE      11              /* significant part of fraction */
 128 #define FRACT_SHIFT     (EXP_SHIFT-FRACT_SIZE)
 129 #define EXP_ADDR(val)   (((val)&EXP_MASK)>>EXP_SHIFT)
 130 #define FRACT_ADDR(val) (((val)&(FRACT_MASK|EXP_LSB))>>FRACT_SHIFT)
 131
 132 #define PR_VEC(a)       a[XX],a[YY],a[ZZ]
 133
 134 #ifdef GMX_SOFTWARE_INVSQRT
 135 extern const unsigned int *  gmx_invsqrt_exptab;
 136 extern const unsigned int *  gmx_invsqrt_fracttab;
 137 #endif
 138
 139
 140 typedef union
 141 {
 142   unsigned int bval;
 143   float fval;
 144 } t_convert;
 145
 146
 147 #ifdef GMX_SOFTWARE_INVSQRT
 148 static real gmx_software_invsqrt(real x)
 149 {
 150   const real  half=0.5;
 151   const real  three=3.0;
 152   t_convert   result,bit_pattern;
 153   unsigned int exp,fract;
 154   real        lu;
 155   real        y;
 156 #ifdef GMX_DOUBLE
 157   real        y2;
 158 #endif
 159
 160   bit_pattern.fval=x;
 161   exp   = EXP_ADDR(bit_pattern.bval);
 162   fract = FRACT_ADDR(bit_pattern.bval);
 163   result.bval=gmx_invsqrt_exptab[exp] | gmx_invsqrt_fracttab[fract];
 164   lu    = result.fval;
 165
 166   y=(half*lu*(three-((x*lu)*lu)));
 167 #ifdef GMX_DOUBLE
 168   y2=(half*y*(three-((x*y)*y)));
 169
 170   return y2;                    /* 10 Flops */
 171 #else
 172   return y;                     /* 5  Flops */
 173 #endif
 174 }
 175 #define gmx_invsqrt(x) gmx_software_invsqrt(x)
 176 #define INVSQRT_DONE
 177 #endif /* gmx_invsqrt */
 178
 179 #ifdef GMX_POWERPC_SQRT
 180 static real gmx_powerpc_invsqrt(real x)
 181 {
 182   const real  half=0.5;
 183   const real  three=3.0;
 184   t_convert   result,bit_pattern;
 185   unsigned int exp,fract;
 186   real        lu;
 187   real        y;
 188 #ifdef GMX_DOUBLE
 189   real        y2;
 190 #endif
 191
 192   lu = __frsqrte((double)x);
 193
 194   y=(half*lu*(three-((x*lu)*lu)));
 195
 196 #if (GMX_POWERPC_SQRT==2)
 197   /* Extra iteration required */
 198   y=(half*y*(three-((x*y)*y)));
 199 #endif
 200
 201 #ifdef GMX_DOUBLE
 202   y2=(half*y*(three-((x*y)*y)));
 203
 204   return y2;                    /* 10 Flops */
 205 #else
 206   return y;                     /* 5  Flops */
 207 #endif
 208 }
 209 #define gmx_invsqrt(x) gmx_powerpc_invsqrt(x)
 210 #define INVSQRT_DONE
 211 #endif /* powerpc_invsqrt */
 212
 213 #ifndef INVSQRT_DONE
 214 #    ifdef GMX_DOUBLE
 215 #        ifdef HAVE_RSQRT
 216 #            define gmx_invsqrt(x)     rsqrt(x)
 217 #        else
 218 #            define gmx_invsqrt(x)     (1.0/sqrt(x))
 219 #        endif
 220 #    else /* single */
 221 #        ifdef HAVE_RSQRTF
 222 #            define gmx_invsqrt(x)     rsqrtf(x)
 223 #        elif defined HAVE_RSQRT
 224 #            define gmx_invsqrt(x)     rsqrt(x)
 225 #        elif defined HAVE_SQRTF
 226 #            define gmx_invsqrt(x)     (1.0/sqrtf(x))
 227 #        else
 228 #            define gmx_invsqrt(x)     (1.0/sqrt(x))
 229 #        endif
 230 #    endif
 231 #endif
 232
 233
 234 static real sqr(real x)
 235 {
 236   return (x*x);
 237 }
 238
 239 static gmx_inline double dsqr(double x)
 240 {
 241   return (x*x);
 242 }
 243
 244 /* Maclaurin series for sinh(x)/x, useful for NH chains and MTTK pressure control
 245    Here, we compute it to 10th order, which might be overkill, 8th is probably enough,
 246    but it's not very much more expensive. */
 247
 248 static gmx_inline real series_sinhx(real x)
 249 {
 250   real x2 = x*x;
 251   return (1 + (x2/6.0)*(1 + (x2/20.0)*(1 + (x2/42.0)*(1 + (x2/72.0)*(1 + (x2/110.0))))));
 252 }
 253
 254 void vecinvsqrt(real in[],real out[],int n);
 255 /* Perform out[i]=1.0/sqrt(in[i]) for n elements */
 256
 257
 258 void vecrecip(real in[],real out[],int n);
 259 /* Perform out[i]=1.0/(in[i]) for n elements */
 260
 261 /* Note: If you need a fast version of vecinvsqrt
 262  * and/or vecrecip, call detectcpu() and run the SSE/3DNow/SSE2/Altivec
 263  * versions if your hardware supports it.
 264  *
 265  * To use those routines, your memory HAS TO BE CACHE-ALIGNED.
 266  * Use snew_aligned(ptr,size,32) to allocate and sfree_aligned to free.
 267  */
 268
 269
 270 static gmx_inline void rvec_add(const rvec a,const rvec b,rvec c)
 271 {
 272   real x,y,z;
 273
 274   x=a[XX]+b[XX];
 275   y=a[YY]+b[YY];
 276   z=a[ZZ]+b[ZZ];
 277
 278   c[XX]=x;
 279   c[YY]=y;
 280   c[ZZ]=z;
 281 }
 282
 283 static gmx_inline void dvec_add(const dvec a,const dvec b,dvec c)
 284 {
 285   double x,y,z;
 286
 287   x=a[XX]+b[XX];
 288   y=a[YY]+b[YY];
 289   z=a[ZZ]+b[ZZ];
 290
 291   c[XX]=x;
 292   c[YY]=y;
 293   c[ZZ]=z;
 294 }
 295
 296 static gmx_inline void ivec_add(const ivec a,const ivec b,ivec c)
 297 {
 298   int x,y,z;
 299
 300   x=a[XX]+b[XX];
 301   y=a[YY]+b[YY];
 302   z=a[ZZ]+b[ZZ];
 303
 304   c[XX]=x;
 305   c[YY]=y;
 306   c[ZZ]=z;
 307 }
 308
 309 static gmx_inline void rvec_inc(rvec a,const rvec b)
 310 {
 311   real x,y,z;
 312
 313   x=a[XX]+b[XX];
 314   y=a[YY]+b[YY];
 315   z=a[ZZ]+b[ZZ];
 316
 317   a[XX]=x;
 318   a[YY]=y;
 319   a[ZZ]=z;
 320 }
 321
 322 static gmx_inline void dvec_inc(dvec a,const dvec b)
 323 {
 324   double x,y,z;
 325
 326   x=a[XX]+b[XX];
 327   y=a[YY]+b[YY];
 328   z=a[ZZ]+b[ZZ];
 329
 330   a[XX]=x;
 331   a[YY]=y;
 332   a[ZZ]=z;
 333 }
 334
 335 static gmx_inline void rvec_sub(const rvec a,const rvec b,rvec c)
 336 {
 337   real x,y,z;
 338
 339   x=a[XX]-b[XX];
 340   y=a[YY]-b[YY];
 341   z=a[ZZ]-b[ZZ];
 342
 343   c[XX]=x;
 344   c[YY]=y;
 345   c[ZZ]=z;
 346 }
 347
 348 static gmx_inline void dvec_sub(const dvec a,const dvec b,dvec c)
 349 {
 350   double x,y,z;
 351
 352   x=a[XX]-b[XX];
 353   y=a[YY]-b[YY];
 354   z=a[ZZ]-b[ZZ];
 355
 356   c[XX]=x;
 357   c[YY]=y;
 358   c[ZZ]=z;
 359 }
 360
 361 static gmx_inline void rvec_dec(rvec a,const rvec b)
 362 {
 363   real x,y,z;
 364
 365   x=a[XX]-b[XX];
 366   y=a[YY]-b[YY];
 367   z=a[ZZ]-b[ZZ];
 368
 369   a[XX]=x;
 370   a[YY]=y;
 371   a[ZZ]=z;
 372 }
 373
 374 static gmx_inline void copy_rvec(const rvec a,rvec b)
 375 {
 376   b[XX]=a[XX];
 377   b[YY]=a[YY];
 378   b[ZZ]=a[ZZ];
 379 }
 380
 381 static gmx_inline void copy_rvecn(rvec *a,rvec *b,int startn, int endn)
 382 {
 383   int i;
 384   for (i=startn;i<endn;i++) {
 385     b[i][XX]=a[i][XX];
 386     b[i][YY]=a[i][YY];
 387     b[i][ZZ]=a[i][ZZ];
 388   }
 389 }
 390
 391 static gmx_inline void copy_dvec(const dvec a,dvec b)
 392 {
 393   b[XX]=a[XX];
 394   b[YY]=a[YY];
 395   b[ZZ]=a[ZZ];
 396 }
 397
 398 static gmx_inline void copy_ivec(const ivec a,ivec b)
 399 {
 400   b[XX]=a[XX];
 401   b[YY]=a[YY];
 402   b[ZZ]=a[ZZ];
 403 }
 404
 405 static gmx_inline void ivec_sub(const ivec a,const ivec b,ivec c)
 406 {
 407   int x,y,z;
 408
 409   x=a[XX]-b[XX];
 410   y=a[YY]-b[YY];
 411   z=a[ZZ]-b[ZZ];
 412
 413   c[XX]=x;
 414   c[YY]=y;
 415   c[ZZ]=z;
 416 }
 417
 418 static gmx_inline void copy_mat(matrix a,matrix b)
 419 {
 420   copy_rvec(a[XX],b[XX]);
 421   copy_rvec(a[YY],b[YY]);
 422   copy_rvec(a[ZZ],b[ZZ]);
 423 }
 424
 425 static gmx_inline void svmul(real a,const rvec v1,rvec v2)
 426 {
 427   v2[XX]=a*v1[XX];
 428   v2[YY]=a*v1[YY];
 429   v2[ZZ]=a*v1[ZZ];
 430 }
 431
 432 static gmx_inline void dsvmul(double a,const dvec v1,dvec v2)
 433 {
 434   v2[XX]=a*v1[XX];
 435   v2[YY]=a*v1[YY];
 436   v2[ZZ]=a*v1[ZZ];
 437 }
 438
 439 static gmx_inline real distance2(const rvec v1,const rvec v2)
 440 {
 441   return sqr(v2[XX]-v1[XX]) + sqr(v2[YY]-v1[YY]) + sqr(v2[ZZ]-v1[ZZ]);
 442 }
 443
 444 static gmx_inline void clear_rvec(rvec a)
 445 {
 446   /* The ibm compiler has problems with inlining this
 447    * when we use a const real variable
 448    */
 449   a[XX]=0.0;
 450   a[YY]=0.0;
 451   a[ZZ]=0.0;
 452 }
 453
 454 static gmx_inline void clear_dvec(dvec a)
 455 {
 456   /* The ibm compiler has problems with inlining this
 457    * when we use a const real variable
 458    */
 459   a[XX]=0.0;
 460   a[YY]=0.0;
 461   a[ZZ]=0.0;
 462 }
 463
 464 static gmx_inline void clear_ivec(ivec a)
 465 {
 466   a[XX]=0;
 467   a[YY]=0;
 468   a[ZZ]=0;
 469 }
 470
 471 static gmx_inline void clear_rvecs(int n,rvec v[])
 472 {
 473 /*  memset(v[0],0,DIM*n*sizeof(v[0][0])); */
 474   int i;
 475
 476   for(i=0; (i<n); i++)
 477     clear_rvec(v[i]);
 478 }
 479
 480 static gmx_inline void clear_mat(matrix a)
 481 {
 482 /*  memset(a[0],0,DIM*DIM*sizeof(a[0][0])); */
 483
 484   const real nul=0.0;
 485
 486   a[XX][XX]=a[XX][YY]=a[XX][ZZ]=nul;
 487   a[YY][XX]=a[YY][YY]=a[YY][ZZ]=nul;
 488   a[ZZ][XX]=a[ZZ][YY]=a[ZZ][ZZ]=nul;
 489 }
 490
 491 static gmx_inline real iprod(const rvec a,const rvec b)
 492 {
 493   return (a[XX]*b[XX]+a[YY]*b[YY]+a[ZZ]*b[ZZ]);
 494 }
 495
 496 static gmx_inline double diprod(const dvec a,const dvec b)
 497 {
 498   return (a[XX]*b[XX]+a[YY]*b[YY]+a[ZZ]*b[ZZ]);
 499 }
 500
 501 static gmx_inline int iiprod(const ivec a,const ivec b)
 502 {
 503   return (a[XX]*b[XX]+a[YY]*b[YY]+a[ZZ]*b[ZZ]);
 504 }
 505
 506 static gmx_inline real norm2(const rvec a)
 507 {
 508   return a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ];
 509 }
 510
 511 static gmx_inline double dnorm2(const dvec a)
 512 {
 513   return a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ];
 514 }
 515
 516 /* WARNING:
 517  * As dnorm() uses sqrt() (which is slow) _only_ use it if you are sure you
 518  * don't need 1/dnorm(), otherwise use dnorm2()*dinvnorm(). */
 519 static gmx_inline double dnorm(const dvec a)
 520 {
 521   return sqrt(diprod(a, a));
 522 }
 523
 524 /* WARNING:
 525  * As norm() uses sqrtf() (which is slow) _only_ use it if you are sure you
 526  * don't need 1/norm(), otherwise use norm2()*invnorm(). */
 527 static gmx_inline real norm(const rvec a)
 528 {
 529   /* This is ugly, but we deliberately do not define gmx_sqrt() and handle the
 530    * float/double case here instead to avoid gmx_sqrt() being accidentally used. */
 531 #ifdef GMX_DOUBLE
 532   return dnorm(a);
 533 #elif defined HAVE_SQRTF
 534   return sqrtf(iprod(a, a));
 535 #else
 536   return sqrt(iprod(a, a));
 537 #endif
 538 }
 539
 540 static gmx_inline real invnorm(const rvec a)
 541 {
 542     return gmx_invsqrt(norm2(a));
 543 }
 544
 545 static gmx_inline real dinvnorm(const dvec a)
 546 {
 547     return gmx_invsqrt(dnorm2(a));
 548 }
 549
 550 /* WARNING:
 551  * Do _not_ use these routines to calculate the angle between two vectors
 552  * as acos(cos_angle(u,v)). While it might seem obvious, the acos function
 553  * is very flat close to -1 and 1, which will lead to accuracy-loss.
 554  * Instead, use the new gmx_angle() function directly.
 555  */
 556 static gmx_inline real
 557 cos_angle(const rvec a,const rvec b)
 558 {
 559   /*
 560    *                  ax*bx + ay*by + az*bz
 561    * cos-vec (a,b) =  ---------------------
 562    *                      ||a|| * ||b||
 563    */
 564   real   cosval;
 565   int    m;
 566   double aa,bb,ip,ipa,ipb,ipab; /* For accuracy these must be double! */
 567
 568   ip=ipa=ipb=0.0;
 569   for(m=0; (m<DIM); m++) {              /* 18           */
 570     aa   = a[m];
 571     bb   = b[m];
 572     ip  += aa*bb;
 573     ipa += aa*aa;
 574     ipb += bb*bb;
 575   }
 576   ipab = ipa*ipb;
 577   if (ipab > 0)
 578     cosval = ip*gmx_invsqrt(ipab);              /*  7           */
 579   else
 580     cosval = 1;
 581                                         /* 25 TOTAL     */
 582   if (cosval > 1.0)
 583     return  1.0;
 584   if (cosval <-1.0)
 585     return -1.0;
 586
 587   return cosval;
 588 }
 589
 590 /* WARNING:
 591  * Do _not_ use these routines to calculate the angle between two vectors
 592  * as acos(cos_angle(u,v)). While it might seem obvious, the acos function
 593  * is very flat close to -1 and 1, which will lead to accuracy-loss.
 594  * Instead, use the new gmx_angle() function directly.
 595  */
 596 static gmx_inline real
 597 cos_angle_no_table(const rvec a,const rvec b)
 598 {
 599   /* This version does not need the invsqrt lookup table */
 600   real   cosval;
 601   int    m;
 602   double aa,bb,ip,ipa,ipb; /* For accuracy these must be double! */
 603
 604   ip=ipa=ipb=0.0;
 605   for(m=0; (m<DIM); m++) {              /* 18           */
 606     aa   = a[m];
 607     bb   = b[m];
 608     ip  += aa*bb;
 609     ipa += aa*aa;
 610     ipb += bb*bb;
 611   }
 612   cosval=ip/sqrt(ipa*ipb);              /* 12           */
 613                                         /* 30 TOTAL     */
 614   if (cosval > 1.0)
 615     return  1.0;
 616   if (cosval <-1.0)
 617     return -1.0;
 618
 619   return cosval;
 620 }
 621
 622
 623 static gmx_inline void cprod(const rvec a,const rvec b,rvec c)
 624 {
 625   c[XX]=a[YY]*b[ZZ]-a[ZZ]*b[YY];
 626   c[YY]=a[ZZ]*b[XX]-a[XX]*b[ZZ];
 627   c[ZZ]=a[XX]*b[YY]-a[YY]*b[XX];
 628 }
 629
 630 static gmx_inline void dcprod(const dvec a,const dvec b,dvec c)
 631 {
 632   c[XX]=a[YY]*b[ZZ]-a[ZZ]*b[YY];
 633   c[YY]=a[ZZ]*b[XX]-a[XX]*b[ZZ];
 634   c[ZZ]=a[XX]*b[YY]-a[YY]*b[XX];
 635 }
 636
 637 /* This routine calculates the angle between a & b without any loss of accuracy close to 0/PI.
 638  * If you only need cos(theta), use the cos_angle() routines to save a few cycles.
 639  * This routine is faster than it might appear, since atan2 is accelerated on many CPUs (e.g. x86).
 640  */
 641 static gmx_inline real
 642 gmx_angle(const rvec a, const rvec b)
 643 {
 644     rvec w;
 645     real wlen,s;
 646
 647     cprod(a,b,w);
 648
 649     wlen  = norm(w);
 650     s     = iprod(a,b);
 651
 652     return atan2(wlen,s);
 653 }
 654
 655 static gmx_inline void mmul_ur0(matrix a,matrix b,matrix dest)
 656 {
 657   dest[XX][XX]=a[XX][XX]*b[XX][XX];
 658   dest[XX][YY]=0.0;
 659   dest[XX][ZZ]=0.0;
 660   dest[YY][XX]=a[YY][XX]*b[XX][XX]+a[YY][YY]*b[YY][XX];
 661   dest[YY][YY]=                    a[YY][YY]*b[YY][YY];
 662   dest[YY][ZZ]=0.0;
 663   dest[ZZ][XX]=a[ZZ][XX]*b[XX][XX]+a[ZZ][YY]*b[YY][XX]+a[ZZ][ZZ]*b[ZZ][XX];
 664   dest[ZZ][YY]=                    a[ZZ][YY]*b[YY][YY]+a[ZZ][ZZ]*b[ZZ][YY];
 665   dest[ZZ][ZZ]=                                        a[ZZ][ZZ]*b[ZZ][ZZ];
 666 }
 667
 668 static gmx_inline void mmul(matrix a,matrix b,matrix dest)
 669 {
 670   dest[XX][XX]=a[XX][XX]*b[XX][XX]+a[XX][YY]*b[YY][XX]+a[XX][ZZ]*b[ZZ][XX];
 671   dest[YY][XX]=a[YY][XX]*b[XX][XX]+a[YY][YY]*b[YY][XX]+a[YY][ZZ]*b[ZZ][XX];
 672   dest[ZZ][XX]=a[ZZ][XX]*b[XX][XX]+a[ZZ][YY]*b[YY][XX]+a[ZZ][ZZ]*b[ZZ][XX];
 673   dest[XX][YY]=a[XX][XX]*b[XX][YY]+a[XX][YY]*b[YY][YY]+a[XX][ZZ]*b[ZZ][YY];
 674   dest[YY][YY]=a[YY][XX]*b[XX][YY]+a[YY][YY]*b[YY][YY]+a[YY][ZZ]*b[ZZ][YY];
 675   dest[ZZ][YY]=a[ZZ][XX]*b[XX][YY]+a[ZZ][YY]*b[YY][YY]+a[ZZ][ZZ]*b[ZZ][YY];
 676   dest[XX][ZZ]=a[XX][XX]*b[XX][ZZ]+a[XX][YY]*b[YY][ZZ]+a[XX][ZZ]*b[ZZ][ZZ];
 677   dest[YY][ZZ]=a[YY][XX]*b[XX][ZZ]+a[YY][YY]*b[YY][ZZ]+a[YY][ZZ]*b[ZZ][ZZ];
 678   dest[ZZ][ZZ]=a[ZZ][XX]*b[XX][ZZ]+a[ZZ][YY]*b[YY][ZZ]+a[ZZ][ZZ]*b[ZZ][ZZ];
 679 }
 680
 681 static gmx_inline void transpose(matrix src,matrix dest)
 682 {
 683   dest[XX][XX]=src[XX][XX];
 684   dest[YY][XX]=src[XX][YY];
 685   dest[ZZ][XX]=src[XX][ZZ];
 686   dest[XX][YY]=src[YY][XX];
 687   dest[YY][YY]=src[YY][YY];
 688   dest[ZZ][YY]=src[YY][ZZ];
 689   dest[XX][ZZ]=src[ZZ][XX];
 690   dest[YY][ZZ]=src[ZZ][YY];
 691   dest[ZZ][ZZ]=src[ZZ][ZZ];
 692 }
 693
 694 static gmx_inline void tmmul(matrix a,matrix b,matrix dest)
 695 {
 696   /* Computes dest=mmul(transpose(a),b,dest) - used in do_pr_pcoupl */
 697   dest[XX][XX]=a[XX][XX]*b[XX][XX]+a[YY][XX]*b[YY][XX]+a[ZZ][XX]*b[ZZ][XX];
 698   dest[XX][YY]=a[XX][XX]*b[XX][YY]+a[YY][XX]*b[YY][YY]+a[ZZ][XX]*b[ZZ][YY];
 699   dest[XX][ZZ]=a[XX][XX]*b[XX][ZZ]+a[YY][XX]*b[YY][ZZ]+a[ZZ][XX]*b[ZZ][ZZ];
 700   dest[YY][XX]=a[XX][YY]*b[XX][XX]+a[YY][YY]*b[YY][XX]+a[ZZ][YY]*b[ZZ][XX];
 701   dest[YY][YY]=a[XX][YY]*b[XX][YY]+a[YY][YY]*b[YY][YY]+a[ZZ][YY]*b[ZZ][YY];
 702   dest[YY][ZZ]=a[XX][YY]*b[XX][ZZ]+a[YY][YY]*b[YY][ZZ]+a[ZZ][YY]*b[ZZ][ZZ];
 703   dest[ZZ][XX]=a[XX][ZZ]*b[XX][XX]+a[YY][ZZ]*b[YY][XX]+a[ZZ][ZZ]*b[ZZ][XX];
 704   dest[ZZ][YY]=a[XX][ZZ]*b[XX][YY]+a[YY][ZZ]*b[YY][YY]+a[ZZ][ZZ]*b[ZZ][YY];
 705   dest[ZZ][ZZ]=a[XX][ZZ]*b[XX][ZZ]+a[YY][ZZ]*b[YY][ZZ]+a[ZZ][ZZ]*b[ZZ][ZZ];
 706 }
 707
 708 static gmx_inline void mtmul(matrix a,matrix b,matrix dest)
 709 {
 710   /* Computes dest=mmul(a,transpose(b),dest) - used in do_pr_pcoupl */
 711   dest[XX][XX]=a[XX][XX]*b[XX][XX]+a[XX][YY]*b[XX][YY]+a[XX][ZZ]*b[XX][ZZ];
 712   dest[XX][YY]=a[XX][XX]*b[YY][XX]+a[XX][YY]*b[YY][YY]+a[XX][ZZ]*b[YY][ZZ];
 713   dest[XX][ZZ]=a[XX][XX]*b[ZZ][XX]+a[XX][YY]*b[ZZ][YY]+a[XX][ZZ]*b[ZZ][ZZ];
 714   dest[YY][XX]=a[YY][XX]*b[XX][XX]+a[YY][YY]*b[XX][YY]+a[YY][ZZ]*b[XX][ZZ];
 715   dest[YY][YY]=a[YY][XX]*b[YY][XX]+a[YY][YY]*b[YY][YY]+a[YY][ZZ]*b[YY][ZZ];
 716   dest[YY][ZZ]=a[YY][XX]*b[ZZ][XX]+a[YY][YY]*b[ZZ][YY]+a[YY][ZZ]*b[ZZ][ZZ];
 717   dest[ZZ][XX]=a[ZZ][XX]*b[XX][XX]+a[ZZ][YY]*b[XX][YY]+a[ZZ][ZZ]*b[XX][ZZ];
 718   dest[ZZ][YY]=a[ZZ][XX]*b[YY][XX]+a[ZZ][YY]*b[YY][YY]+a[ZZ][ZZ]*b[YY][ZZ];
 719   dest[ZZ][ZZ]=a[ZZ][XX]*b[ZZ][XX]+a[ZZ][YY]*b[ZZ][YY]+a[ZZ][ZZ]*b[ZZ][ZZ];
 720 }
 721
 722 static gmx_inline real det(matrix a)
 723 {
 724   return ( a[XX][XX]*(a[YY][YY]*a[ZZ][ZZ]-a[ZZ][YY]*a[YY][ZZ])
 725           -a[YY][XX]*(a[XX][YY]*a[ZZ][ZZ]-a[ZZ][YY]*a[XX][ZZ])
 726           +a[ZZ][XX]*(a[XX][YY]*a[YY][ZZ]-a[YY][YY]*a[XX][ZZ]));
 727 }
 728
 729 static gmx_inline void m_add(matrix a,matrix b,matrix dest)
 730 {
 731   dest[XX][XX]=a[XX][XX]+b[XX][XX];
 732   dest[XX][YY]=a[XX][YY]+b[XX][YY];
 733   dest[XX][ZZ]=a[XX][ZZ]+b[XX][ZZ];
 734   dest[YY][XX]=a[YY][XX]+b[YY][XX];
 735   dest[YY][YY]=a[YY][YY]+b[YY][YY];
 736   dest[YY][ZZ]=a[YY][ZZ]+b[YY][ZZ];
 737   dest[ZZ][XX]=a[ZZ][XX]+b[ZZ][XX];
 738   dest[ZZ][YY]=a[ZZ][YY]+b[ZZ][YY];
 739   dest[ZZ][ZZ]=a[ZZ][ZZ]+b[ZZ][ZZ];
 740 }
 741
 742 static gmx_inline void m_sub(matrix a,matrix b,matrix dest)
 743 {
 744   dest[XX][XX]=a[XX][XX]-b[XX][XX];
 745   dest[XX][YY]=a[XX][YY]-b[XX][YY];
 746   dest[XX][ZZ]=a[XX][ZZ]-b[XX][ZZ];
 747   dest[YY][XX]=a[YY][XX]-b[YY][XX];
 748   dest[YY][YY]=a[YY][YY]-b[YY][YY];
 749   dest[YY][ZZ]=a[YY][ZZ]-b[YY][ZZ];
 750   dest[ZZ][XX]=a[ZZ][XX]-b[ZZ][XX];
 751   dest[ZZ][YY]=a[ZZ][YY]-b[ZZ][YY];
 752   dest[ZZ][ZZ]=a[ZZ][ZZ]-b[ZZ][ZZ];
 753 }
 754
 755 static gmx_inline void msmul(matrix m1,real r1,matrix dest)
 756 {
 757   dest[XX][XX]=r1*m1[XX][XX];
 758   dest[XX][YY]=r1*m1[XX][YY];
 759   dest[XX][ZZ]=r1*m1[XX][ZZ];
 760   dest[YY][XX]=r1*m1[YY][XX];
 761   dest[YY][YY]=r1*m1[YY][YY];
 762   dest[YY][ZZ]=r1*m1[YY][ZZ];
 763   dest[ZZ][XX]=r1*m1[ZZ][XX];
 764   dest[ZZ][YY]=r1*m1[ZZ][YY];
 765   dest[ZZ][ZZ]=r1*m1[ZZ][ZZ];
 766 }
 767
 768 static gmx_inline void m_inv_ur0(matrix src,matrix dest)
 769 {
 770   double tmp = src[XX][XX]*src[YY][YY]*src[ZZ][ZZ];
 771   if (fabs(tmp) <= 100*GMX_REAL_MIN)
 772     gmx_fatal(FARGS,"Can not invert matrix, determinant is zero");
 773
 774   dest[XX][XX] = 1/src[XX][XX];
 775   dest[YY][YY] = 1/src[YY][YY];
 776   dest[ZZ][ZZ] = 1/src[ZZ][ZZ];
 777   dest[ZZ][XX] = (src[YY][XX]*src[ZZ][YY]*dest[YY][YY]
 778                   - src[ZZ][XX])*dest[XX][XX]*dest[ZZ][ZZ];
 779   dest[YY][XX] = -src[YY][XX]*dest[XX][XX]*dest[YY][YY];
 780   dest[ZZ][YY] = -src[ZZ][YY]*dest[YY][YY]*dest[ZZ][ZZ];
 781   dest[XX][YY] = 0.0;
 782   dest[XX][ZZ] = 0.0;
 783   dest[YY][ZZ] = 0.0;
 784 }
 785
 786 static gmx_inline void m_inv(matrix src,matrix dest)
 787 {
 788   const real smallreal = (real)1.0e-24;
 789   const real largereal = (real)1.0e24;
 790   real  deter,c,fc;
 791
 792   deter = det(src);
 793   c     = (real)1.0/deter;
 794   fc    = (real)fabs(c);
 795
 796   if ((fc <= smallreal) || (fc >= largereal))
 797     gmx_fatal(FARGS,"Can not invert matrix, determinant = %e",deter);
 798
 799   dest[XX][XX]= c*(src[YY][YY]*src[ZZ][ZZ]-src[ZZ][YY]*src[YY][ZZ]);
 800   dest[XX][YY]=-c*(src[XX][YY]*src[ZZ][ZZ]-src[ZZ][YY]*src[XX][ZZ]);
 801   dest[XX][ZZ]= c*(src[XX][YY]*src[YY][ZZ]-src[YY][YY]*src[XX][ZZ]);
 802   dest[YY][XX]=-c*(src[YY][XX]*src[ZZ][ZZ]-src[ZZ][XX]*src[YY][ZZ]);
 803   dest[YY][YY]= c*(src[XX][XX]*src[ZZ][ZZ]-src[ZZ][XX]*src[XX][ZZ]);
 804   dest[YY][ZZ]=-c*(src[XX][XX]*src[YY][ZZ]-src[YY][XX]*src[XX][ZZ]);
 805   dest[ZZ][XX]= c*(src[YY][XX]*src[ZZ][YY]-src[ZZ][XX]*src[YY][YY]);
 806   dest[ZZ][YY]=-c*(src[XX][XX]*src[ZZ][YY]-src[ZZ][XX]*src[XX][YY]);
 807   dest[ZZ][ZZ]= c*(src[XX][XX]*src[YY][YY]-src[YY][XX]*src[XX][YY]);
 808 }
 809
 810 static gmx_inline void mvmul(matrix a,const rvec src,rvec dest)
 811 {
 812   dest[XX]=a[XX][XX]*src[XX]+a[XX][YY]*src[YY]+a[XX][ZZ]*src[ZZ];
 813   dest[YY]=a[YY][XX]*src[XX]+a[YY][YY]*src[YY]+a[YY][ZZ]*src[ZZ];
 814   dest[ZZ]=a[ZZ][XX]*src[XX]+a[ZZ][YY]*src[YY]+a[ZZ][ZZ]*src[ZZ];
 815 }
 816
 817 static gmx_inline void mvmul_ur0(matrix a,const rvec src,rvec dest)
 818 {
 819   dest[ZZ]=a[ZZ][XX]*src[XX]+a[ZZ][YY]*src[YY]+a[ZZ][ZZ]*src[ZZ];
 820   dest[YY]=a[YY][XX]*src[XX]+a[YY][YY]*src[YY];
 821   dest[XX]=a[XX][XX]*src[XX];
 822 }
 823
 824 static gmx_inline void tmvmul_ur0(matrix a,const rvec src,rvec dest)
 825 {
 826   dest[XX]=a[XX][XX]*src[XX]+a[YY][XX]*src[YY]+a[ZZ][XX]*src[ZZ];
 827   dest[YY]=                  a[YY][YY]*src[YY]+a[ZZ][YY]*src[ZZ];
 828   dest[ZZ]=                                    a[ZZ][ZZ]*src[ZZ];
 829 }
 830
 831 static gmx_inline void unitv(const rvec src,rvec dest)
 832 {
 833   real linv;
 834
 835   linv=gmx_invsqrt(norm2(src));
 836   dest[XX]=linv*src[XX];
 837   dest[YY]=linv*src[YY];
 838   dest[ZZ]=linv*src[ZZ];
 839 }
 840
 841 static gmx_inline void unitv_no_table(const rvec src,rvec dest)
 842 {
 843   real linv;
 844
 845   linv=1.0/sqrt(norm2(src));
 846   dest[XX]=linv*src[XX];
 847   dest[YY]=linv*src[YY];
 848   dest[ZZ]=linv*src[ZZ];
 849 }
 850
 851 static void calc_lll(rvec box,rvec lll)
 852 {
 853   lll[XX] = 2.0*M_PI/box[XX];
 854   lll[YY] = 2.0*M_PI/box[YY];
 855   lll[ZZ] = 2.0*M_PI/box[ZZ];
 856 }
 857
 858 static gmx_inline real trace(matrix m)
 859 {
 860   return (m[XX][XX]+m[YY][YY]+m[ZZ][ZZ]);
 861 }
 862
 863 static gmx_inline real _divide_err(real a,real b,const char *file,int line)
 864 {
 865     if (fabs(b) <= GMX_REAL_MIN)
 866         gmx_fatal(FARGS,"Dividing by zero, file %s, line %d",file,line);
 867     return a/b;
 868 }
 869
 870 static gmx_inline int _mod(int a,int b,char *file,int line)
 871 {
 872   if(b==0)
 873     gmx_fatal(FARGS,"Modulo zero, file %s, line %d",file,line);
 874   return a % b;
 875 }
 876
 877 /* Operations on multidimensional rvecs, used e.g. in edsam.c */
 878 static void m_rveccopy(int dim, rvec *a, rvec *b)
 879 {
 880     /* b = a */
 881     int i;
 882
 883     for (i=0; i<dim; i++)
 884         copy_rvec(a[i],b[i]);
 885 }
 886
 887 /*computer matrix vectors from base vectors and angles */
 888 static void matrix_convert(matrix box, rvec vec, rvec angle)
 889 {
 890     svmul(DEG2RAD,angle,angle);
 891     box[XX][XX] = vec[XX];
 892     box[YY][XX] = vec[YY]*cos(angle[ZZ]);
 893     box[YY][YY] = vec[YY]*sin(angle[ZZ]);
 894     box[ZZ][XX] = vec[ZZ]*cos(angle[YY]);
 895     box[ZZ][YY] = vec[ZZ]
 896                          *(cos(angle[XX])-cos(angle[YY])*cos(angle[ZZ]))/sin(angle[ZZ]);
 897     box[ZZ][ZZ] = sqrt(sqr(vec[ZZ])
 898                        -box[ZZ][XX]*box[ZZ][XX]-box[ZZ][YY]*box[ZZ][YY]);
 899 }
 900
 901 #define divide_err(a,b) _divide_err((a),(b),__FILE__,__LINE__)
 902 #define mod(a,b)    _mod((a),(b),__FILE__,__LINE__)
 903
 904 #ifdef __cplusplus
 905 }
 906 #endif
 907
 908
 909 #endif  /* _vec_h */