/*
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"
/*
extern void
CopyBlochPad(struct Grid *outre, struct Grid *outim,
struct Grid *inre, struct Grid *inim,
struct Dvec kvec, struct Ivec pad);
extern void
CopyfromPadded(struct Grid *output, struct Grid *input, struct Ivec pad);
*/
void apply_I(const WL_ComplexColumn &input, WL_ComplexColumn &output)
{
#ifdef DFT_PROFILING
timerOn(34); // Turn on I timer
counterIncr(12); // increment I counter
#endif // DFT_PROFILING
if (input.representation == REALSPACE)
dft_log("**> Trying to take I of realspace column.\n");
// Set the representation and quantum numbers of the output.
output.representation=REALSPACE;
// Declare workspace variables for operators.
real *workspace1, *workspace2;
// Allocate workspaces
workspace1= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace1", "apply_I");
workspace2= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace2", "apply_I");
//switch between algorithms for sparse/dense
if(input.embedding==output.embedding){
Gridspec *gridspec;
// choose the appropriate <sparse|dense> gridspec
switch(input.embedding){
case SPARSE:
gridspec = input.basis_spec->gridspec;
break;
case DENSE:
gridspec = input.basis_spec->gridspec_dense;
break;
default:
die("apply_I: Unknown input embedding %d\n",input.embedding);
}
Grid *temp = mkgridnew(gridspec,NULL);
// I also works in place.
GetRe(temp,input); zeroghostrec(temp);
I(temp,workspace1,workspace2,input.basis->lwork);
PutRe(temp,output);
GetIm(temp,input); zeroghostrec(temp);
I(temp,workspace1,workspace2,input.basis->lwork);
PutIm(temp,output);
killgridnew(temp);
} else if (input.embedding==SPARSE && output.embedding==DENSE){
// die("I(%s, %s) not implemented yet\n",
// input.getembedding(),output.getembedding());
Gridspec *ingridspec=input.basis_spec->gridspec;
Gridspec *outgridspec=input.basis_spec->gridspec_dense;
Grid *colre = mkgridnew(ingridspec,NULL);
Grid *colim = mkgridnew(ingridspec,NULL);
Grid *Icolre = mkgridnew(outgridspec,NULL);
Grid *Icolim = mkgridnew(outgridspec,NULL);
GetRe(colre, input); zeroghostrec(colre);
GetIm(colim, input); zeroghostrec(colim);
struct Dvec kvec;
if(input.basis->qnum){ // hack it chargedensities don't have qnum
kvec.x=input.basis->qnum->kvec.v[0];
kvec.y=input.basis->qnum->kvec.v[1];
kvec.z=input.basis->qnum->kvec.v[2];
} else kvec=NULLDVEC;
I_S2D(Icolre,Icolim,colre, colim, kvec,\
workspace1,workspace2,input.basis->lwork);
PutRe(Icolre, output);
PutIm(Icolim, output);
killgridnew(colre); killgridnew(colim);
killgridnew(Icolre); killgridnew(Icolim);
} else
die("Don't know how to I(%s, %s)\n",
input.getembedding(), output.getembedding());
myfree(workspace1); myfree(workspace2);
#ifdef DFT_PROFILING
timerOff(34); // Turn off I timer
#endif // DFT_PROFILING
}
void
apply_J(const WL_ComplexColumn &input, WL_ComplexColumn &output)
{
#ifdef DFT_PROFILING
timerOn(35); // Turn on J timer
counterIncr(14); // J counter
#endif // DFT_PROFILING
if (input.representation == COEFFSPACE)
dft_log("**> Trying to take J of coeffspace column.\n");
// Set the representation and quantum numbers of the output.
output.representation=COEFFSPACE;
if(output.embedding!=input.embedding)
die("Cannot do apply_J(%s, %s)! Write me!\n", input.getembedding(), output.getembedding());
// output.kvec=input.kvec;
// Declare workspace variables for operators.
real *workspace1, *workspace2;
// Allocate workspaces
workspace1= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace1", "apply_J");
workspace2= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace2", "apply_J");
if(input.embedding!=output.embedding)
die("Embeddings must match for J(col, col)\n");
Gridspec *gridspec;
// choose the appropriate gridspec
switch(input.embedding){
case SPARSE:
gridspec = input.basis_spec->gridspec;
break;
case DENSE:
gridspec = input.basis_spec->gridspec_dense;
break;
default:
die("apply_J: Unknown input embedding %d\n",input.embedding);
}
Grid *temp= mkgridnew(gridspec,NULL);
// J appears to work in-place on the input.
GetRe(temp,input);
adjustdown(temp);
J(temp,workspace1,workspace2,input.basis->lwork);
PutRe(temp,output);
GetIm(temp,input);
adjustdown(temp);
J(temp,workspace1,workspace2,input.basis->lwork);
PutIm(temp,output);
killgridnew(temp);
myfree(workspace1); myfree(workspace2);
#ifdef DFT_PROFILING
timerOff(35); // Turn off J timer
#endif // DFT_PROFILING
}
void
apply_Idag(const WL_ComplexColumn &input, WL_ComplexColumn &output)
{
#ifdef DFT_PROFILING
timerOn(36); // Turn on Idag timer
counterIncr(13); // Idag counter
#endif // DFT_PROFILING
if (input.representation==COEFFSPACE)
dft_log("**> Trying to take Idag of coeffspace column.\n");
// Set the representation and quantum numbers of the output.
output.representation=COEFFSPACE;
// Declare workspace variables for operators.
real *workspace1, *workspace2;
// Allocate workspaces
workspace1= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace1", "apply_Idag");
workspace2= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace2", "apply_Idag");
//switch between algorithms for sparse/dense
if(input.embedding==output.embedding){
Gridspec *gridspec;
// choose the appropriate <sparse|dense> gridspec
switch(input.embedding){
case SPARSE:
gridspec = input.basis_spec->gridspec;
break;
case DENSE:
gridspec = input.basis_spec->gridspec_dense;
break;
default:
die("apply_Idag: Unknown input embedding %d\n",input.embedding);
}
Grid *temp = mkgridnew(gridspec,NULL);
// Idag works in-place.
GetRe(temp,input);
adjustdown(temp); // shouldn't be necessary
Idag(temp,workspace1,workspace2,input.basis->lwork);
PutRe(temp,output);
GetIm(temp,input);
adjustdown(temp); // shouldn't be necessary
Idag(temp,workspace1,workspace2,input.basis->lwork);
PutIm(temp,output);
// release the temp memory
killgridnew(temp);
}else if (input.embedding==DENSE && output.embedding==SPARSE){
// die("apply_Idag(%s, %s) not implemented yet\n",
// input.getembedding(), output.getembedding());
Gridspec *ingridspec=input.basis_spec->gridspec_dense;
Gridspec *outgridspec=input.basis_spec->gridspec;
Grid *colre = mkgridnew(ingridspec,NULL);
Grid *colim = mkgridnew(ingridspec,NULL);
Grid *Icolre = mkgridnew(outgridspec,NULL);
Grid *Icolim = mkgridnew(outgridspec,NULL);
GetRe(colre, input);
GetIm(colim, input);
struct Dvec kvec;
if(input.basis->qnum){ // hack it chargedensities don't have qnum
kvec.x=input.basis->qnum->kvec.v[0];
kvec.y=input.basis->qnum->kvec.v[1];
kvec.z=input.basis->qnum->kvec.v[2];
} else kvec=NULLDVEC;
Idag_D2S(Icolre,Icolim,colre, colim, kvec,\
workspace1,workspace2,input.basis->lwork);
PutRe(Icolre, output);
PutIm(Icolim, output);
killgridnew(colre); killgridnew(colim);
killgridnew(Icolre); killgridnew(Icolim);
} else
die("Don't know how to Idag(%s, %s)\n",
input.getembedding(), output.getembedding());
myfree(workspace1); myfree(workspace2);
#ifdef DFT_PROFILING
timerOff(36); // Turn off Idag timer
#endif // DFT_PROFILING
}
void
apply_Jdag(const WL_ComplexColumn &input, WL_ComplexColumn &output)
{
#ifdef DFT_PROFILING
timerOn(37); // Turn on Jdag timer
counterIncr(15); // Jdag counter
#endif // DFT_PROFILING
if (input.representation==REALSPACE)
dft_log("**> Trying to take Jdag of realspace column.\n");
// Set the representation and quantum numbers of the output.
output.representation=REALSPACE;
// Allocate workspaces
real *workspace1= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace1", "apply_Jdag");
real *workspace2= (real *) mymalloc(input.basis->lwork*sizeof(real),
"workspace2", "apply_Jdag");
if(input.embedding!=output.embedding)
die("Embeddings must match for Jdag(col, col)\n");
Gridspec *gridspec;
// choose the appropriate gridspec
switch(input.embedding){
case SPARSE:
gridspec = input.basis_spec->gridspec;
break;
case DENSE:
gridspec = input.basis_spec->gridspec_dense;
break;
default:
die("apply_Jdag: Unknown input embedding %d\n",input.embedding);
}
Grid *temp = mkgridnew(gridspec,NULL);
// Jdag works in-place.
GetRe(temp,input);
zeroghostrec(temp); //shouldn't be necessary
Jdag(temp,workspace1,workspace2,input.basis->lwork);
PutRe(temp,output);
GetIm(temp,input);
zeroghostrec(temp); //shouldn't be necessary
Jdag(temp,workspace1,workspace2,input.basis->lwork);
PutIm(temp,output);
// release the temp memory
killgridnew(temp);
myfree(workspace1); myfree(workspace2);
#ifdef DFT_PROFILING
timerOff(37); // Turn off Jdag timer
#endif // DFT_PROFILING
}
void
apply_D(int d,int dag,const WL_ComplexColumn &in, WL_ComplexColumn &out)
{
// out=in;
// return;
/* Sanity checks */
if(dag<0 || dag>1)
die("Illegal value of input variable dag in apply_D()!\n");
if(d<0 || dag>2)
die("Illegal value of input variable d in apply_D()!\n");
/* Process real/coeff space issues */
if (dag==0)
if (in.representation == COEFFSPACE)
printf("Trying to take Di of coeffspace column.\n");
else
out.representation=REALSPACE;
else
if (in.representation == COEFFSPACE)
printf("Trying to take Didag of coeffspace column.\n");
else
out.representation=REALSPACE;
if(out.embedding!=in.embedding)
die("Cannot do apply_D(%s, %s)! Write me!\n", in.getembedding(), out.getembedding());
// Declare workspace variables for operators.
real *work1, *work2;
// Allocate workspaces
work1= (real *) mymalloc(in.basis->lwork*sizeof(real),
"work1", "apply_D");
work2= (real *) mymalloc(in.basis->lwork*sizeof(real),
"work2", "apply_D");
if(in.embedding!=out.embedding)
die("Embeddings must match for apply_D(col, col)\n");
Gridspec *gridspec;
// choose the appropriate gridspec
switch(in.embedding){
case SPARSE:
gridspec = in.basis_spec->gridspec;
break;
case DENSE:
gridspec = in.basis_spec->gridspec_dense;
break;
default:
die("apply_D: Unknown input embedding %d\n",in.embedding);
}
/* This seems like a LOT of extra grids. Could improve this by
reprogramming D and Ddag to work in place */
Grid *temp1= mkgridnew(gridspec,NULL);
Grid *temp2= mkgridnew(gridspec,NULL);
Grid *work= mkgridnew(gridspec,NULL);
GetRe(temp1,in);
if (dag==0) {
J(temp1,work1,work2,in.basis->lwork);
D(d,temp2,temp1,work,work1,work2,in.basis->lwork);
}
else {
Ddag(d,temp2,temp1,work,work1,work2,in.basis->lwork);
Jdag(temp2,work1,work2,in.basis->lwork);
}
PutRe(temp2,out);
GetIm(temp1,in);
if (dag==0) {
J(temp1,work1,work2,in.basis->lwork);
D(d,temp2,temp1,work,work1,work2,in.basis->lwork);
}
else {
Ddag(d,temp2,temp1,work,work1,work2,in.basis->lwork);
Jdag(temp2,work1,work2,in.basis->lwork);
}
PutIm(temp2,out);
killgridnew(temp1);
killgridnew(temp2);
killgridnew(work);
myfree(work1); myfree(work2);
}
syntax highlighted by Code2HTML, v. 0.9.1