/* 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. */ /* Actual Implementation of the Transforms I,J, Idag,Jdag for Plane Wave Basis They all act on a single column */ #include "header.h" // We do not have the operations on single argument, because there are // issues with those, on what do they have to return //I takes us from Fourier to Real space void apply_I(const PW_ComplexColumn &col, PW_ComplexColumn &Icol) { #ifdef DFT_PROFILING timerOn(34); // Turn on I timer counterIncr(12); // increment I counter #endif // DFT_PROFILING if(col.representation == REALSPACE) die("Trying to call I on realspace column! What r u thinking?!\n"); // call a map in case we need it. if not, it just does a copy col.map(Icol); FFT3D(Icol.basis->basis_spec->Nx, Icol.basis->basis_spec->Ny, Icol.basis->basis_spec->Nz, Icol.data.d , 1); /* Get rid of FFT scaling and divide by sqrt(Vol) */ real s = (real)1.0/sqrt(Icol.basis->basis_spec->lattice->unit_cell_volume); Icol *= s; //keep track of representation flag Icol.representation=REALSPACE; #ifdef DFT_PROFILING timerOff(34); // Turn off I timer #endif // DFT_PROFILING } /* The J operator on column_bundles: this takes from r to G space */ void apply_J(const PW_ComplexColumn &col, PW_ComplexColumn &Jcol) { #ifdef DFT_PROFILING timerOn(35); // Turn on J timer counterIncr(14); // J counter #endif // DFT_PROFILING if(col.representation == COEFFSPACE) die("Trying to call J on coeffspace column! What r u thinking?!\n"); //ipd: put some extra code, to avoid this copy if both are full size columns PW_ComplexColumn Jcolfull(col); FFT3D(Jcolfull.basis->basis_spec->Nx, Jcolfull.basis->basis_spec->Ny, Jcolfull.basis->basis_spec->Nz, Jcolfull.data.d , -1); Jcolfull.map(Jcol); /* Get rid of FFT scaling and divide by sqrt(Vol) */ real s = sqrt(Jcol.basis->basis_spec->lattice->unit_cell_volume); Jcol *= s; Jcol.representation = COEFFSPACE; #ifdef DFT_PROFILING timerOff(35); // Turn off J timer #endif // DFT_PROFILING } /* applies Idag=(Nx*Ny*Nz)/Vol*J to Y and puts it into IdagY */ void apply_Idag(const PW_ComplexColumn &col, PW_ComplexColumn &Idagcol) { #ifdef DFT_PROFILING timerOn(36); // Turn on Idag timer counterIncr(13); // Idag counter #endif // DFT_PROFILING apply_J(col, Idagcol); real s = (real)col.basis->basis_spec->NxNyNz/col.basis->basis_spec->lattice->unit_cell_volume; Idagcol *= s; #ifdef DFT_PROFILING timerOff(36); // Turn off Idag timer #endif // DFT_PROFILING } /* We use Jdag = Vol/(Nx*Ny*Nz)*I */ void apply_Jdag(const PW_ComplexColumn &col, PW_ComplexColumn &Jdagcol) { #ifdef DFT_PROFILING timerOn(37); // Turn on Jdag timer counterIncr(15); // Jdag counter #endif // DFT_PROFILING if(col.representation == REALSPACE) die("Trying to call Jdag on realspace column! What r u thinking?!\n"); apply_I(col, Jdagcol); real s = col.basis->basis_spec->lattice->unit_cell_volume/(real)col.basis->basis_spec->NxNyNz; Jdagcol *= s; #ifdef DFT_PROFILING timerOff(37); // Turn off Jdag timer #endif // DFT_PROFILING } // // This is a specialized routine that is used when applying the Hamiltonian // (apply_Hsp() in PHLO.c) which saves one vector worth of memory // by destroying the contents of Y. // /* Idag = Nx*Ny*Nz*FFT(-)/sqrt(Vol) = (Nx*Ny*Nz)*J/Vol. * This routine applies Idag to Y and places the result into IdagY, but * in the process destroys the contents of Y (i.e. it does FFT in place * on Y). */ // we probably can't do in-place transforms on the top level /* template */ /* void */ /* apply_Idag_inplace(PW_ComplexColumn *col, PW_ComplexColumn *Idagcol) */ /* { */ /* dft_log("apply_Idag inplace: calling complex FFTW3D\n"); */ /* } */ /* The Dd operator on column_bundles (maps r to r space): d=x,y,z dag=0,1 (0 for Dd, 1 for Dd-dagger) */ void apply_D(int d,int dag,const PW_ComplexColumn &in, PW_ComplexColumn &out) { /* Sanity checks */ if(in.representation == COEFFSPACE) die("Trying to call Dx on coeffspace column! What r u thinking?!\n"); if ( in.length!=out.length ) die("apply_D(in,out) called with different sizes for in and out\n"); /* Extract pointers to relevant info */ Basis *basis = in.basis; Lattice *lattice = basis->basis_spec->lattice; if (basis == NULL) die("apply_D(in,out) called with in.basis == NULL\n"); /* Copy data immediately into Dxcol and work on it there... */ out=in; /* A true FFT inverse */ int Nx,Ny,Nz,Nx2,Ny2,Nz2,NyNz,NxNyNz; 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; NxNyNz=basis->basis_spec->NxNyNz; FFT3D(Nx, Ny, Nz, out.data.d, -1); /* Apply gradient in Fourier space */ double Gd,Gd0,Gd1,Gd2; Gd0 = lattice->G.m[0][d]; Gd1 = lattice->G.m[1][d]; Gd2 = lattice->G.m[2][d]; /* Convert (0,1) flag to proper sign factor for loop */ if (dag<0 || dag>1) die("apply_D called with inappropriate value for flag dag\n"); dag=1-2*dag; /* Loop over Fourier space, zeroing Niquist edges!!! */ for (int i=-Nx2; i < Nx2; i++) for (int j=-Ny2; j < Ny2; j++) for (int k=-Nz2; k < Nz2; k++) { int 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; Gd = i*Gd0+j*Gd1+k*Gd2; if (i==-Nx2 || j==-Ny2 || k==-Nz2 ) out.data.d[index] = complex(0.,0.); else out.data.d[index] *= complex(0.,dag*Gd); } /* Forward FFT (no need to scale -here-!) */ FFT3D(Nx, Ny, Nz, out.data.d, 1); out.representation=in.representation; }