/* 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 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 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); }