/* DFT++ is a density functional package developed by the research group of Professor Tomas Arias Copyright 1996-2003 Sohrab Ismail-Beigi This file is part of DFT++. DFT++ is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. DFT++ is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with DFT++; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Please see the file CREDITS for a list of authors. For academic users, we request that publications using results obtained with this software reference "New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000). and, if using the wavelet basis, further reference "Multiresolution analysis of electronic structure: semicardinal and wavelet bases," T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999). and "Robust ab initio calculation of condensed matter: transparent convergence through semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003). For your convenience, preprints of the above articles may be obtained from http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively. */ /* * Set of 3-vector and 3x3-matrix classes and routines written by * Gabor Csanyi, Mar 1997. * */ /* $Id: lin3.h,v 1.5.2.3 2003/05/29 18:54:27 ivan Exp $ */ #ifndef DFT_LIN3_H #define DFT_LIN3_H // "real" must already been defined, say as double. #include #include /* * 3-vector class: * * Operators: +,-,+=,-=,= are standard. * * between two vectors is the dot-product. * * between a vector and a scalar is scaling. * / between a vector and a scalar is scaling. * ^ between two vectors is the cross-product. * */ class vector3 { /* variables */ public: real v[3]; /* constructors */ vector3() { v[0] = v[1] = v[2] = 0.0; } vector3(real aa, real bb, real cc) { v[0] = aa; v[1] = bb; v[2] = cc; } vector3(const vector3 &vec) /* copy constructor */ {v[0] = vec.v[0]; v[1] = vec.v[1]; v[2] = vec.v[2]; } /* various operators: */ /* assignment */ vector3 operator=(const vector3 &vec){ int i; for(i=0;i<3;i++) v[i] = vec.v[i]; return(*this); } /* assign to constant */ vector3 operator=(const real val) { for (int i=0; i<3; i++) v[i] = val; return(*this); } /* unary - */ vector3 operator-() const { vector3 vec( -v[0], -v[1], -v[2] ); return vec; } /* sum two vectors */ vector3 operator+(const vector3 &vec) const { vector3 nvec( v[0]+vec.v[0], v[1]+vec.v[1], v[2]+vec.v[2] ); return(nvec); } /* accumulate sum of two vectors */ vector3 operator+=(const vector3 &vec){ int i; for(i=0;i<3;i++) v[i] += vec.v[i]; return(*this); } /* subtract vectors */ vector3 operator-(const vector3 &vec) const { vector3 nvec( v[0]-vec.v[0], v[1]-vec.v[1], v[2]-vec.v[2] ); return(nvec); } /* do -= */ vector3 operator-=(const vector3 &vec){ int i; for(i=0;i<3;i++) v[i] -= vec.v[i]; return(*this); } /* 3-vector dot product */ real operator*(const vector3 &vec) const { return(v[0]*vec.v[0]+v[1]*vec.v[1]+v[2]*vec.v[2]); } /* 3-vector cross product */ vector3 operator^(const vector3 &vec) const { vector3 nvec( v[1]*vec.v[2] - v[2]*vec.v[1], v[2]*vec.v[0] - v[0]*vec.v[2], v[0]*vec.v[1] - v[1]*vec.v[0] ); return(nvec); } /* scaling vector by real number (various versions): s*v */ friend vector3 operator*(real s, const vector3 &vec); /* v*s */ vector3 operator*(real s) const { vector3 vec( s*v[0], s*v[1], s*v[2] ); return(vec); } /* v *= s */ vector3 operator*=(real s) { for (int i=0; i<3; i++) v[i] *= s; return(*this); } /* v/s */ vector3 operator/(real s) const { vector3 vec( v[0]/s, v[1]/s, v[2]/s ); return(vec); } /* v/w , matlab style: v./w */ vector3 operator/(vector3 w) const { vector3 vec( v[0]/w.v[0], v[1]/w.v[1], v[2]/w.v[2] ); return(vec); } /* print 3-vector */ void print(Output *out, char *format){ int i; out->printf("[ "); for (i=0; i < 3; i++) out->printf(format,v[i]); out->printf(" ]\n"); } /* test equality */ int operator==(vector3 w) const { for(int i=0; i < 3; i++) if(v[i] != w.v[i]) return 0; return 1; } }; /* s*vec */ inline vector3 operator*(real s, const vector3 &vec) { vector3 nvec( s*vec.v[0], s*vec.v[1], s*vec.v[2] ); return(nvec); } /* Cartesian length of 3 vector */ inline real abs(const vector3 &vec) { return(sqrt(vec.v[0]*vec.v[0]+vec.v[1]*vec.v[1]+vec.v[2]*vec.v[2])); } /* Cartesian length^2 of 3 vector */ inline real abs2(const vector3 &vec) { return(vec.v[0]*vec.v[0]+vec.v[1]*vec.v[1]+vec.v[2]*vec.v[2]); } /* * 3D integer vector class: * ipd: i need this for the WL symemtries - just a structure * */ /* #define DFT_DVEC2IVEC_TOL .1 class ivector3 { public: int x, y, z; ivector3(){ x=y=z=0; } ivector3(class vector3 vec){ if(fabs(fabs(x=(int)vec.v[0])-fabs(vec.v[0]))>DFT_DVEC2IVEC_TOL) die("you don't want to assign %lf to %d\n", vec.v[0], x); if(fabs(fabs(y=(int)vec.v[1])-fabs(vec.v[1]))>DFT_DVEC2IVEC_TOL) die("you don't want to assign %lf to %d\n", vec.v[1], y); if(fabs(fabs(z=(int)vec.v[2])-fabs(vec.v[2]))>DFT_DVEC2IVEC_TOL) die("you don't want to assign %lf to %d\n", vec.v[2], z); } }; */ /* * 3x3 matrix class. * * Operators: +,+=,-,-=,= are standard * *,/ with real number scales the matrix * *,*= matrix matrix multiply * [i] returns i'th column of matrix (vector3) * ~ returns transpose of matrix */ class matrix3 { /* data */ public: real m[3][3]; /* constructors:*/ /* init to zero */ matrix3() { int i, j; for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j]=0.0; } /* init to diagonal */ matrix3(real d1, real d2, real d3) { m[0][0] = d1; m[1][1] = d2, m[2][2] = d3; m[0][1] = m[0][2] = m[1][0] = m[1][2] = m[2][0] = m[2][1] = 0.0; } /* init row by row */ matrix3(real m11, real m12, real m13, real m21, real m22, real m23, real m31, real m32, real m33){ m[0][0] = m11; m[0][1] = m12; m[0][2] = m13; m[1][0] = m21; m[1][1] = m22; m[1][2] = m23; m[2][0] = m31; m[2][1] = m32; m[2][2] = m33; } /* join together 3 column vectors */ matrix3(vector3 v1, vector3 v2, vector3 v3) { m[0][0] = v1.v[0]; m[0][1] = v2.v[0]; m[0][2] = v3.v[0]; m[1][0] = v1.v[1]; m[1][1] = v2.v[1]; m[1][2] = v3.v[1]; m[2][0] = v1.v[2]; m[2][1] = v2.v[2]; m[2][2] = v3.v[2]; } /* copy constructor */ matrix3(const matrix3 &n){ int i, j; for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j]=n.m[i][j]; } /* operators: */ /* assignment */ matrix3 operator=(const matrix3 &n){ int i, j; for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j]=n.m[i][j]; return(*this); } /* comparison */ int operator==(const matrix3 &n) const { int i, j; for (i=0; i<3; i++) { for (j=0; j<3; j++) { if (m[i][j] != n.m[i][j]) { return 0; } } } return 1; } /* ~comparison */ int operator!=(const matrix3 &n) const { int i, j; for (i=0; i<3; i++) { for (j=0; j<3; j++) { if (m[i][j] != n.m[i][j]) { return 1; } } } return 0; } /* unary - */ matrix3 operator-() const { matrix3 n; int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)n.m[i][j]=-m[i][j]; return(n); } /* add two matrices */ matrix3 operator+(const matrix3 &n) const { matrix3 t; int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)t.m[i][j]=n.m[i][j]+m[i][j]; return(t); } /* += a matrix */ matrix3 operator+=(const matrix3 &n){ int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j]+=n.m[i][j]; return(*this); } /* subtract two matrices */ matrix3 operator-(const matrix3 &n) const { matrix3 t; int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)t.m[i][j]=m[i][j]-n.m[i][j]; return(t); } /* -= */ matrix3 operator-=(const matrix3 &n){ int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j]-=n.m[i][j]; return(*this); } /* matrix multiply */ matrix3 operator*(const matrix3 &n) const { matrix3 t; int i,j, k; for(i=0;i<3;i++)for(j=0;j<3;j++){ t.m[i][j] = 0.0; for(k = 0; k < 3; k++) t.m[i][j] += m[i][k]*n.m[k][j]; } return(t); } /* matrix times vector */ vector3 operator*(const vector3 &vec) const { vector3 nvec( m[0][0]*vec.v[0]+m[0][1]*vec.v[1]+m[0][2]*vec.v[2], m[1][0]*vec.v[0]+m[1][1]*vec.v[1]+m[1][2]*vec.v[2], m[2][0]*vec.v[0]+m[2][1]*vec.v[1]+m[2][2]*vec.v[2] ); return(nvec); } /* scale matrix by real number (various versions): m*s */ matrix3 operator*(real s) const { matrix3 t; int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)t.m[i][j]=s*m[i][j]; return(t); } /* m *= s */ matrix3 operator*=(real s){ int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j] *= s; return(*this); } /* s*m */ friend matrix3 operator*(real s,const matrix3 &n); /* matrix mult. *= */ matrix3 operator*=(const matrix3 &n) { matrix3 t; int i,j, k; for(i=0;i<3;i++)for(j=0;j<3;j++){ t.m[i][j] = 0.0; for(k = 0; k < 2; k++) t.m[i][j] += n.m[i][k]*m[k][j]; } for(i=0;i<3;i++)for(j=0;j<3;j++)m[i][j]=t.m[i][j]; return(*this); } /* scale by 1/s real number */ matrix3 operator/(real s) const { matrix3 n; int i, j; for(i=0;i<3;i++)for(j=0;j<3;j++)n.m[i][j]=m[i][j]/s; return(n); } /* extract i'th column of a matrix */ vector3 operator[](int i) const { vector3 vec; vec.v[0] = m[0][i]; vec.v[1] = m[1][i]; vec.v[2] = m[2][i]; return(vec); } /* transpose */ matrix3 operator~() const { matrix3 n; int i,j; for(i = 0; i < 3; i++) for(j = 0; j < 3; j++) n.m[i][j] = m[j][i]; return(n); } /* print it */ void print(Output *out,char *format) { int i,j; for (i=0; i < 3; i++) { out->printf("[ "); for (j=0; j < 3; j++) out->printf(format,m[i][j]); out->printf(" ]\n"); } } }; /* s*N */ inline matrix3 operator*(real s,const matrix3 &n) { matrix3 t; int i,j; for(i=0;i<3;i++)for(j=0;j<3;j++)t.m[i][j]=s*n.m[i][j]; return(t); } /* Trace */ inline real trace3(const matrix3 &n) { return(n.m[0][0]+n.m[1][1]+n.m[2][2]); } /* Determinant */ inline real det3(const matrix3 &n) { return(n.m[0][0]*(n.m[1][1]*n.m[2][2]-n.m[1][2]*n.m[2][1])+ n.m[0][1]*(n.m[1][2]*n.m[2][0]-n.m[1][0]*n.m[2][2])+ n.m[0][2]*(n.m[1][0]*n.m[2][1]-n.m[1][1]*n.m[2][0])); } /* Inverse */ inline matrix3 inv3(const matrix3 &n) { matrix3 nt; nt = ~n; matrix3 minv(nt[1]^nt[2], nt[2]^nt[0], nt[0]^nt[1]); minv *= 1.0/det3(n); return(minv); } #endif // DFT_LIN3_H