MolLib/Types.h

   1 /*
   2  *  Types.h
   3  *
   4  *  Created by Toshi Nagata on 06/03/11.
   5  *  Copyright 2006-2008 Toshi Nagata. All rights reserved.
   6  *
   7  This program is free software; you can redistribute it and/or modify
   8  it under the terms of the GNU General Public License as published by
   9  the Free Software Foundation version 2 of the License.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15  */
  16
  17 #ifndef __Types_h__
  18 #define __Types_h__
  19
  20 #include <stdlib.h>
  21 #include <stdio.h>
  22 #include <math.h>
  23
  24 #if defined(__WXMAC__) || defined(__CMDMAC__)
  25 /*  On Mac OS X, CLAPACK is in Accelerate.framework  */
  26 #include <vecLib/cblas.h>
  27 #include <vecLib/clapack.h>
  28 #else
  29 #include <f2c.h>
  30 #include <blaswrap.h>
  31 #include <clapack.h>
  32 #endif
  33
  34 /*  Get the eigenvalue/eigenvector for a real symmetric matrix (3x3)  */
  35 #if !defined(__WXMAC__) && !defined(__CMDMAC__)
  36 typedef integer        __CLPK_integer;
  37 typedef logical        __CLPK_logical;
  38 typedef real           __CLPK_real;
  39 typedef doublereal     __CLPK_doublereal;
  40 #endif
  41
  42 #ifdef __cplusplus
  43 extern "C" {
  44 #endif
  45
  46 #define STUB extern
  47
  48 #ifndef kRad2Deg
  49 #define kRad2Deg  (180./3.14159265358979)
  50 #endif
  51 #ifndef kDeg2Rad
  52 #define kDeg2Rad  (3.14159265358979 / 180.)
  53 #endif
  54 #ifndef kLog10
  55 #define kLog10 2.3025851
  56 #endif
  57 #ifndef kArcCot15
  58 #define kArcCot15 3.73205075
  59 #endif
  60
  61 #undef kBohr2Angstrom
  62 #define kBohr2Angstrom 0.529177249
  63
  64 #undef kAngstrom2Bohr
  65 #define kAngstrom2Bohr 1.8897259885
  66
  67 typedef double Double;
  68 typedef int Int;
  69 typedef unsigned char Byte;
  70 typedef unsigned int UInt;
  71
  72 #define kInvalidIndex -99999999  /*  Used for terminating integer array  */
  73
  74 #define kInvalidFloat 1e30
  75
  76 typedef struct Vector { Double x, y, z; } Vector;
  77 typedef struct Quat { Double x, y, z, w; } Quat;   /*  A quaternion  */
  78 typedef Double Mat33[9];  /*  Columns first!  */
  79
  80 /*  Mat33 (rot) and Vector (translation) */
  81 /*  Transform[0..8] are the same as Mat33, and Transform[9..11] are the fourth column. So that, the
  82     matrix elements are in the following order;
  83     {a11, a21, a31, a12, a22, a32, a13, a23, a33, a14, a24, a34}  */
  84 /*  (Changed from row-first Mat33 on 2011.12.7.)  */
  85 typedef Double Transform[12];
  86 typedef Double *TransformPtr;
  87 /*typedef Double Matrix[4][4]; */
  88
  89 #define VecAdd(v3, v1, v2) ((v3).x=(v1).x+(v2).x, (v3).y=(v1).y+(v2).y, (v3).z=(v1).z+(v2).z)
  90 #define VecInc(v1, v2) ((v1).x+=(v2).x, (v1).y+=(v2).y, (v1).z+=(v2).z)
  91 #define VecSub(v3, v1, v2) ((v3).x=(v1).x-(v2).x, (v3).y=(v1).y-(v2).y, (v3).z=(v1).z-(v2).z)
  92 #define VecDec(v1, v2) ((v1).x-=(v2).x, (v1).y-=(v2).y, (v1).z-=(v2).z)
  93 #define VecScale(v2, v1, r) ((v2).x=(v1).x*(r), (v2).y=(v1).y*(r), (v2).z=(v1).z*(r))
  94 #define VecScaleSelf(v1, r) ((v1).x*=(r), (v1).y*=(r), (v1).z*=(r))
  95 #define VecScaleInc(v2, v1, r) ((v2).x+=(v1).x*(r), (v2).y+=(v1).y*(r), (v2).z+=(v1).z*(r))
  96 #define VecLength2(v1) ((v1).x*(v1).x+(v1).y*(v1).y+(v1).z*(v1).z)
  97 #define VecLength(v1) (sqrt(VecLength2(v1)))
  98 #define VecDot(v1, v2) ((v1).x*(v2).x+(v1).y*(v2).y+(v1).z*(v2).z)
  99 #define VecCross(v3, v1, v2) ((v3).x=(v1).y*(v2).z-(v1).z*(v2).y, (v3).y=(v1).z*(v2).x-(v1).x*(v2).z, (v3).z=(v1).x*(v2).y-(v1).y*(v2).x)
 100 #define VecZero(v1) ((v1).x=(v1).y=(v1).z=0.0)
 101 #define VecIndex(vp, i) (((Double *)(vp))[i])
 102
 103 /*  Vector utility functions  */
 104 int NormalizeVec(Vector *vdst, const Vector *vsrc);  /*  Returns non-zero when zero-vector */
 105
 106 void MatrixRotation(Mat33 dst, const Vector *axis, Double angle);
 107 void MatrixVec(Vector *dst, const Mat33 mat, const Vector *src);
 108 void MatrixMul(Mat33 dst, const Mat33 src1, const Mat33 src2);
 109 void MatrixTranspose(Mat33 dst, const Mat33 src);
 110 void MatrixScale(Mat33 dst, const Mat33 src, Double factor);
 111 Double MatrixDeterminant(const Mat33 src);
 112 int MatrixInvert(Mat33 dst, const Mat33 src);  /*  Return non-zero when determinant is zero */
 113 int MatrixSymDiagonalize(Mat33 mat, Double *out_values, Vector *out_vectors);
 114 void MatrixGeneralRotation(Mat33 dst, const Vector *v1, const Vector *v2, const Vector *v3);
 115
 116 void TransformVec(Vector *dst, const Transform tf, const Vector *src);
 117 void TransformMul(Transform dst, const Transform src1, const Transform src2);
 118 Double TransformDeterminant(const Transform src);
 119 int TransformInvert(Transform dst, const Transform src);  /*  Return non-zero when determinant is zero */
 120 void TransformTranspose(Transform dst, const Transform src);
 121
 122 void TransformForInversion(Transform dst, const Vector *center);
 123 int  TransformForReflection(Transform dst, const Vector *axis, const Vector *center);
 124 int  TransformForRotation(Transform dst, const Vector *axis, Double angle, const Vector *center);
 125
 126 extern Transform gIdentityTransform;
 127
 128 /*  Wrapper struct for CLAPACK routines  */
 129 typedef struct LAMatrix {
 130         __CLPK_integer row, column;
 131         __CLPK_doublereal data[1];
 132 } LAMatrix;
 133
 134 LAMatrix *LAMatrixAllocTempMatrix(int row, int column);
 135 void LAMatrixReleaseTempMatrix(LAMatrix *mat);
 136
 137 LAMatrix *LAMatrixNew(int row, int column);
 138 void LAMatrixRelease(LAMatrix *mat);
 139 LAMatrix *LAMatrixResize(LAMatrix *mat, int row, int column);
 140 LAMatrix *LAMatrixNewFromMatrix(const LAMatrix *mat);
 141 /*  mat3 = scale1 * mat1 * mat2 + scale2 * mat3  */
 142 /*  If trans1/trans2 is non-zero, mat1/mat2 is transposed before multiplication  */
 143 void LAMatrixMul(int trans1, int trans2, double scale1, const LAMatrix *mat1, const LAMatrix *mat2, double scale2, LAMatrix *mat3);
 144 int LAMatrixInvert(LAMatrix *mat1, const LAMatrix *mat2);
 145 Double LAMatrixDeterminant(const LAMatrix *mat);
 146 void LAMatrixTranspose(LAMatrix *mat1, const LAMatrix *mat2);
 147 int LAMatrixSymDiagonalize(LAMatrix *vec, LAMatrix *mat1, const LAMatrix *mat2);
 148 int LAMatrixSingularValueDecomposition(LAMatrix *matU, LAMatrix *matW, LAMatrix *matV, const LAMatrix *mat);
 149
 150 /*  Utility functions  */
 151 void SetPanicFunc(void (*func)(const char *, ...));
 152 void SetWarningFunc(void (*func)(const char *, ...));
 153 extern void (*gPanicFunc)(const char *, ...);
 154 extern void (*gWarningFunc)(const char *, ...);
 155 #define Panic (*gPanicFunc)
 156 #define Warning (*gWarningFunc)
 157
 158 void PanicByOutOfMemory(const char *msg);
 159 void PanicByInternalError(const char *msg, const char *file, int line);
 160
 161 #define MALLOC_CHECK(p, s) ((p) == NULL ? PanicByOutOfMemory(s) : (void)0)
 162 #define NULL_CHECK(p, s) ((p) == NULL ? PanicByInternalError((s), __FILE__, __LINE__) : (void)0)
 163
 164 void *AssignArray(void *base, Int *count, int item_size, int idx, const void *value);
 165 void *NewArray(void *base, Int *count, int item_size, int nitems);
 166 void *InsertArray(void *base, Int *count, int item_size, int idx, int nitems, const void *value);
 167 void *DeleteArray(void *base, Int *count, int item_size, int idx, int nitems, void *outValue);
 168 int ReadLine(char *buf, int size, FILE *stream, int *lineNumber);
 169 int ReadFormat(const char *str, const char *fmt,...);
 170
 171 #ifdef __cplusplus
 172 }
 173 #endif
 174
 175 #endif /* __Types_h__ */