/*
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.
*/
/*************************************************************************
*
* Defines the PlaneWave operations (mostly for the Solid State case)
*
*************************************************************************/
#include "header.h"
//constructors
PW_ComplexColumn::PW_ComplexColumn()
{
fat = TRUE;
length = 0;
basis_spec = NULL;
basis = NULL;
representation = REALSPACE; //default
embedding = SPHERE;
}
PW_ComplexColumn::PW_ComplexColumn(PW_Basis *b, char s[10])
{
if(!strcmp(s,"fat"))
fat = TRUE;
else if (!strcmp(s,"slim"))
fat = FALSE;
else
die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");
basis = b;
representation = COEFFSPACE; //default
if( b != NULL)
basis_spec = b->basis_spec;
}
PW_ComplexColumn::PW_ComplexColumn(PW_Basis *b, int space, char s[10])
{
if(!strcmp(s,"fat"))
fat = TRUE;
else if (!strcmp(s,"slim"))
fat = FALSE;
else
die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");
basis = b;
representation = space;
if( b != NULL)
basis_spec = b->basis_spec;
}
//copy constructor w/ extra argument
PW_ComplexColumn::PW_ComplexColumn(const PW_ComplexColumn &col,
char s[10])
{
if(!strcmp(s,"fat"))
fat = TRUE;
else if (!strcmp(s,"slim"))
fat = FALSE;
else
die("PW_Column(): Identity crisis - don't know if i'm slim or fat\n");
length = col.length;
basis_spec = col.basis_spec;
basis = col.basis;
representation = col.representation;
embedding = col.embedding;
if(fat) {
data.init(length);
data = col.data;
}
else
die("You can't copy a slim column, fatso!\n");
}
/*
//copy constructor
PW_ComplexColumn::PW_ComplexColumn(const PW_ComplexColumn &col)
{
fat = col.fat;
length = col.length;
basis_spec = col.basis_spec;
basis = col.basis;
representation = col.representation;
embedding = col.embedding;
if(fat) {
data.init(length);
data = col.data;
}
else
die("You can't copy a slim column, fatso!\n");
}
*/
// destructor
PW_ComplexColumn::~PW_ComplexColumn()
{
if(fat)
data.free();
}
void
PW_ComplexColumn::init_wavefunction(PW_Basis &b)
{
length = b.nbasis;
basis_spec = b.basis_spec;
basis = &b;
representation = COEFFSPACE;
embedding = SPHERE;
}
void
PW_ComplexColumn::init_scalarfield(PW_Basis &b, int space)
{
basis = &b;
basis_spec = b.basis_spec;
length = basis_spec->NxNyNz;
embedding = BOX;
representation = space;
if(fat)
data.init(length);
else
die("You cannot initialize a slim column to contain a scalar field\n");
}
void
PW_ComplexColumn::print()
{
int j;
for (j = 0; j < length; j++)
dft_log("%f\t%f\n", data.d[j].x, data.d[j].y);
}
char *
PW_ComplexColumn::getrepresentation()
{
static char realspacename[] = "RealSpace";
static char coeffspacename[] = "CoeffSpace";
switch(representation){
case REALSPACE:
return realspacename;
case COEFFSPACE:
return coeffspacename;
default:
die("Unknown representation: %d\n", representation);
}
return NULL;
}
char *
PW_ComplexColumn::getembedding()
{
static char spherename[] = "Sphere";
static char boxname[] = "Box";
switch(embedding){
case SPHERE:
return spherename;
case BOX:
return boxname;
default:
die("Unknown embedding: %d\n", embedding);
}
return NULL;
}
// LINEAR ALGEBRA - just call the array algebra
void
PW_ComplexColumn::operator=(const PW_ComplexColumn &col)
{
if(length != col.length)
die("length != col.length in PW_ComplexColumn::operator=(): %d %d\n", length, col.length);
if(embedding != col.embedding)
die("embedding != col.embedding in PW_ComplexColumn::operator=()\n");
representation = col.representation;
data = col.data;
}
void
PW_ComplexColumn::operator+=(const PW_ComplexColumn &col)
{
if(length != col.length)
die("length != col.length in PW_ComplexColumn::operator+=(): %d %d\n", length, col.length);
if(representation != col.representation)
die("representation != col.repr. in PW_ComplexColumn::operator+=()\n");
if(embedding != col.embedding)
die("embedding != col.embedding in PW_ComplexColumn::operator+=()\n");
data += col.data;
}
void
PW_ComplexColumn::operator*=(const PW_ComplexColumn &col)
{
if(length != col.length)
die("length != col.length in PW_ComplexColumn::operator*=() %d %d\n", length, col.length);
if(representation != col.representation)
die("representation != col.repr. in PW_ComplexColumn::operator*=()\n");
if(embedding != col.embedding)
die("embedding != col.embedding in PW_ComplexColumn::operator*=()\n");
data *= col.data;
}
void
PW_ComplexColumn::operator-=(const PW_ComplexColumn &col)
{
if(length != col.length)
die("length != col.length in PW_ComplexColumn::operator-=()%d %d\n", length, col.length);
if(representation != col.representation)
die("representation != col.repr. in PW_ComplexColumn::operator-=()\n");
if(embedding != col.embedding)
die("embedding != col.embedding in PW_ComplexColumn::operator-=()\n");
data -= col.data;
}
void
PW_ComplexColumn::operator*=(const real r)
{
data*=r;
}
void
PW_ComplexColumn::operator*=(const complex c)
{
data*=c;
}
complex
PW_ComplexColumn::operator^(const PW_ComplexColumn &col)
{
return dot(this->data, col.data);
}
// BASIS SPECIFIC
/* Fill the column with random gaussian distributed numbers
* with zero mean and std. dev. = 1/(1+((0.5*G^2)/0.75)^6), where 0.5*G^2
* is the kinetic energy of the basis function in question.
* The reason for this strange formula is that we want high-energy
* plane-waves to have lower initial weight than low-energy one. The
* cutoff of 0.75 is because following the Jones-scheme business for
* the 2-atom Si cell, we want to keep the (0 0 0), (1 1 1), and (2 0 0)
* vectors (units are 2*pi/a). The energy of the (2 0 0) in Hartrees
* and for a=5.43 Angstorms is 0.7499 (miracle!?!?), and the next
* highest planewave is (2 2 0) with energy 1.4998. So I just chose
* some function which turns over at 0.75 relatively sharply.
*
* -- SIB
*/
void
PW_ComplexColumn::randomize()
{
int j;
real std,KE,t, ttt;
vector3 kplusG;
// ipd: hmm, sth's fishy here -should have it work for anything
// this basically limits us to work on spheres
//if (basis->nbasis != length)
// die("\nPW_ComplexColumn::randomize() nbasis != length!!!\n\n");
switch (representation){
case REALSPACE:
//just put random stuff in there
data.randomize();
break;
case COEFFSPACE:
if(embedding==SPHERE)
for(j=0; j < basis->nbasis; j++){
kplusG.v[0] = (real)basis->Gx[j] + basis->qnum->kvec.v[0];
kplusG.v[1] = (real)basis->Gy[j] + basis->qnum->kvec.v[1];
kplusG.v[2] = (real)basis->Gz[j] + basis->qnum->kvec.v[2];
KE = 0.5*(kplusG*(basis->basis_spec->lattice->GGT*kplusG));
t = KE/0.75;
//std = 1.0/(1.0+t*t*t*t*t*t);
ttt = t*t*t;
t = ttt*ttt;
std = 1.0/(1.0+t);
data.d[j].x = gauss(std);
data.d[j].y = gauss(std);
}
else //embedding == BOX
//i'll put uniform noise for now
data.randomize();
//gabor, looks we need this just for the tests - you can put
//k-dep if you want, i'm leaving to you the resurection ...
/* die("Don't know how to randomize %s %s\n", */
/* getrepresentation(), getembedding()); */
break;
default:
die("unknown representation %s in randomize\n", getrepresentation());
}
}
void
PW_ComplexColumn::setmodones()
{
int ii;
real theta = rand(); // just a random phase
complex one(cos(theta), sin(theta));
for(ii=0; ii<length; ii++)
data.d[ii] = one;
//PW_ComplexColumn out(*this); // create a copy column
if(this->representation == COEFFSPACE)
apply_J(*this, *this);
return;
}
void
PW_ComplexColumn::IdentityMap(PW_ComplexColumn &col) const
{
if(this != &col)
col = *this;
}
void
PW_ComplexColumn::SphereToBoxMap(PW_ComplexColumn &col) const
{
int *index = this->basis->index;
int j;
//ipd: considert col.data.zero_out();
// Fill in part of G-space taken up by the basis with correct values
for (j=0; j < col.basis->basis_spec->NxNyNz; j++){
// dft_log("zeroing %d\n", j);
col.data.d[j] = 0.0;
}
for (j=0; j < basis->nbasis; j++){
//dft_log("mapping %d\n", j);
col.data.d[index[j]] = this->data.d[j];
}
}
void
PW_ComplexColumn::BoxToSphereMap(PW_ComplexColumn &col) const
{
int *index = col.basis->index;
int j;
// Fill in part of G-space taken up by the basis with correct values
for (j=0; j < col.basis->nbasis; j++){
//dft_log("UNmapping %d\n", j);
col.data.d[j] = this->data.d[index[j]];
}
}
//ipd: adding the boxtospheremap and the mapflag functionality
void
PW_ComplexColumn::map(PW_ComplexColumn &mappedCol) const
{
switch(embedding){
case BOX:
switch(mappedCol.embedding){
case BOX:
this->IdentityMap(mappedCol);
break;
case SPHERE:
this->BoxToSphereMap(mappedCol);
break;
default:
die("Invalid embedding in PW_ComplexColumn::map() : %d\n", mappedCol.embedding);
}
break;
case SPHERE:
switch(mappedCol.embedding){
case SPHERE:
this->IdentityMap(mappedCol);
break;
case BOX:
this->SphereToBoxMap(mappedCol);
break;
default:
die("Invalid embedding in PW_ComplexColumn::map() : %d\n", mappedCol.embedding);
}
break;
default:
die("Invalid map in PW_ComplexColumn::map() : %d\n", embedding);
}
}
void
PW_ComplexColumn::read(char *filename)
{
data.read(filename);
}
void
PW_ComplexColumn::write(char *filename)
{
data.write(filename);
}
// adds the square of the absolute value
void
add_scale_abs2(const complex &c, const PW_ComplexColumn & in,
PW_ComplexColumn &out)
{
add_scale_abs2(c, in.data, out.data);
}
//returns the absolute squares of the input column
PW_ComplexColumn abs2(PW_ComplexColumn &in)
{
int i;
PW_ComplexColumn out(in);
for(i=0; i< in.length; i++)
out.data.d[i]=abs2(in.data.d[i]);
return out;
}
PW_ComplexColumn
pointwise_mult(const PW_ComplexColumn &a,
const PW_ComplexColumn &b)
{
if(a.length != b.length)
die("a.length != b.length in pointwise_mult!\n");
PW_ComplexColumn ab(a);
int i;
for(i=0; i<a.length; i++)
ab.data.d[i] *= b.data.d[i];
return ab;
}
// UNUSED CODE
/* Fill the column_bundle with random gaussian distributed numbers
* with zero mean and std. dev. = 1/(1+((0.5*G^2)/0.75)^6), where 0.5*G^2
* is the kinetic energy of the basis function in question.
* The reason for this strange formula is that we want high-energy
* plane-waves to have lower initial weight than low-energy one. The
* cutoff of 0.75 is because following the Jones-scheme business for
* the 2-atom Si cell, we want to keep the (0 0 0), (1 1 1), and (2 0 0)
* vectors (units are 2*pi/a). The energy of the (2 0 0) in Hartrees
* and for a=5.43 Angstorms is 0.7499 (miracle!?!?), and the next
* highest planewave is (2 2 0) with energy 1.4998. So I just chose
* some function which turns over at 0.75 relatively sharply.
*/
/* void */
/* PW_Column<RealArray>::randomize(void) */
/* { */
/* register int i,j; */
/* register real std,KE,t,ttt; */
/* vector3 kplusG; */
/* if (basis->nbasis != length) */
/* die("\nPW_Column<RA>::randomize() nbasis != length!!!\n\n"); */
/* for(j=0; j < basis->nbasis; j++){ */
/* kplusG.v[0] = (real)basis->Gx[j] + qnum->kvec.v[0]; */
/* kplusG.v[1] = (real)basis->Gy[j] + qnum->kvec.v[1]; */
/* kplusG.v[2] = (real)basis->Gz[j] + qnum->kvec.v[2]; */
/* KE = 0.5*(kplusG*(basis->lattice->GGT*kplusG)); */
/* t = KE/0.75; */
/* //std = 1.0/(1.0+t*t*t*t*t*t); */
/* ttt = t*t*t; */
/* t = ttt*ttt; */
/* std = 1.0/(1.0+t); */
/* this->data.d[j] = gauss(std); */
/* } */
/* } */
// this will be different, but let's just leave it like this for now,
// until we figure out the gamma point case
/* void */
/* PW_Column<RealArray>::init_wavefunction(PW_Basis &b) */
/* { */
/* basis = &b; */
/* basis_spec = b.basis_spec; */
/* map = &(PW_Column<RealArray>::SphereToBoxMap); */
/* } */
syntax highlighted by Code2HTML, v. 0.9.1