/*
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.
*/
#include "header.h"
void apply_O(const PW_ComplexColumn &c, PW_ComplexColumn &Oc)
{
#ifdef DFT_PROFILING
timerOn(29); // turn on O timer
#endif // DFT_PROFILING
if(&c!=&Oc)
Oc = c;
#ifdef DFT_PROFILING
timerOff(29); // turn off O timer
#endif // DFT_PROFILING
}
void apply_Obar(const PW_ComplexColumn &c, PW_ComplexColumn &Obarc)
{
if(c.representation == REALSPACE)
die("Trying to take Obar of realspace column!\n");
if(&c != &Obarc)
Obarc = c;
Obarc.data.d[0].x = 0.0;
return;
}
void apply_L(const PW_ComplexColumn &c, PW_ComplexColumn &Lc)
{
#ifdef DFT_PROFILING
timerOn(31); // turn on L timer
#endif // DFT_PROFILING
if(c.representation == REALSPACE)
die("Trying to take L of realspace column!\n");
Lc.representation=c.representation;
int i,j,k, nbasis;
int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz, index;
real kplusGx,kplusGy,kplusGz,kx,ky,kz;
real GGTxx,GGTyy,GGTzz,GGTxy,GGTxz,GGTyz;
real kplusG2, G2;
Basis *basis = c.basis;
Lattice *lattice = basis->basis_spec->lattice;
if (basis == NULL)
die("apply_L() called with c.basis == NULL\n");
if ( c.length!=Lc.length )
die("apply_L() called with different sizes for c and Lc\n");
GGTxx = lattice->GGT.m[0][0];
GGTyy = lattice->GGT.m[1][1];
GGTzz = lattice->GGT.m[2][2];
GGTxy = lattice->GGT.m[0][1];
GGTxz = lattice->GGT.m[0][2];
GGTyz = lattice->GGT.m[1][2];
/*
* we need to decide if we have to compute G vectors, or if they
* are stored for us.
*/
if(basis->Gx == NULL || basis->Gy == NULL || basis->Gz == NULL){
// we compute G vectors
Nx = basis->basis_spec->Nx; Nx2 = Nx/2;
Ny = basis->basis_spec->Ny; Ny2 = Ny/2;
Nz = basis->basis_spec->Nz; Nz2 = Nz/2;
NyNz = Ny*Nz;
for (i=-Nx2; i < Nx2; i++)
for (j=-Ny2; j < Ny2; j++)
for (k=-Nz2; k < Nz2; k++)
{
index = 0;
if (k < 0) index += k+Nz; else index += k;
if (j < 0) index += Nz*(j+Ny); else index += Nz*j;
if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
G2 = i*i*GGTxx +
j*j*GGTyy +
k*k*GGTzz +
2*(i*j*GGTxy +
i*k*GGTxz +
j*k*GGTyz);
Lc.data.d[index] = -G2*c.data.d[index];
}
} else {
// G vectors are stored
// make sure we have a quantum number
if(basis->qnum == NULL)
die("In apply_L() G vectors are not stored and qnum == NULL\n");
kx = basis->qnum->kvec.v[0];
ky = basis->qnum->kvec.v[1];
kz = basis->qnum->kvec.v[2];
nbasis = basis->nbasis;
for (j=0; j < nbasis; j++){
kplusGx = kx + basis->Gx[j];
kplusGy = ky + basis->Gy[j];
kplusGz = kz + basis->Gz[j];
kplusG2 = (kplusGx*kplusGx*GGTxx +
kplusGy*kplusGy*GGTyy +
kplusGz*kplusGz*GGTzz +
2*(kplusGx*kplusGy*GGTxy +
kplusGx*kplusGz*GGTxz +
kplusGy*kplusGz*GGTyz ));
Lc.data.d[j] = -kplusG2*c.data.d[j];
}
}
#ifdef DFT_PROFILING
timerOff(31); // turn off L timer
#endif // DFT_PROFILING
}
/*
* The inverse of the Laplacian on vectors: the is obviously modulo
* G=0 where there are infinities. For planewaves, we have -1/|G|^2.
*/
void apply_invL(const PW_ComplexColumn &c, PW_ComplexColumn &invLc)
{
#ifdef DFT_PROFILING
timerOn(32); // turn on invL timer
#endif // DFT_PROFILING
int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz;
int i,j,k,index;
real G2,GGTxx,GGTxy,GGTxz,GGTyy,GGTyz,GGTzz;
BasisSpec *spec;
Lattice *lattice;
if (c.basis == 0)
die("invL(vector) called with v.basis==0\n");
if (c.length != c.basis_spec->NxNyNz)
die("invL(vector) called with vector.n != FFT box size\n");
if ((c.basis_spec->Nx%2)!=0 || (c.basis_spec->Ny%2)!=0 ||
(c.basis_spec->Nz%2)!=0)
die("invL(vector) called with basis->Nx,Ny,or Nz not a multiple of 2\n");
spec=c.basis_spec;
lattice = spec->lattice;
Nx = spec->Nx; Nx2 = Nx/2;
Ny = spec->Ny; Ny2 = Ny/2;
Nz = spec->Nz; Nz2 = Nz/2;
NyNz = Ny*Nz;
GGTxx = lattice->GGT.m[0][0];
GGTyy = lattice->GGT.m[1][1];
GGTzz = lattice->GGT.m[2][2];
GGTxy = lattice->GGT.m[0][1];
GGTxz = lattice->GGT.m[0][2];
GGTyz = lattice->GGT.m[1][2];
for (i=-Nx2; i < Nx2; i++)
for (j=-Ny2; j < Ny2; j++)
for (k=-Nz2; k < Nz2; k++)
{
/* G = 0 case */
if (i==0 && j==0 && k==0)
invLc.data.d[0] = 0.0;
else
{
index = 0;
if (k < 0) index += k+Nz; else index += k;
if (j < 0) index += Nz*(j+Ny); else index += Nz*j;
if (i < 0) index += NyNz*(i+Nx); else index += NyNz*i;
G2 = i*i*GGTxx +
j*j*GGTyy +
k*k*GGTzz +
2*(i*j*GGTxy +
i*k*GGTxz +
j*k*GGTyz);
invLc.data.d[index] = -c.data.d[index]/G2;
}
}
#ifdef DFT_PROFILING
timerOff(32); // turn off invL timer
#endif // DFT_PROFILING
return;
}
syntax highlighted by Code2HTML, v. 0.9.1