/*
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 <stdio.h>
#include <math.h>
/*
* 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
syntax highlighted by Code2HTML, v. 0.9.1