/*
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<class myArray1, class myArray2> */
/* 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;
}
syntax highlighted by Code2HTML, v. 0.9.1