/*
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.
*/
/*
* locpseudopot.c: Sohrab Ismail-Beigi
* Written Aug. 1996, modified Jan. 97, Feb. 97
*
* Routines that give the atomic pseudopotential as a function of r or k,
* and Vpseudogrid() finds the crystal pseudopotential on a specified grid.
* There are various convergence approximations... see below
* for details.
*
* complex structfact()
*/
#include "header.h"
/*
* S(G) = sum_[j=1,natoms] { exp(-i*2*pi*G.r_j) }: the infamous
* structure factor. Here G is in reciprocal lattice coords and hence
* given by integers; the atomic positions are in lattice coordinates.
*/
complex structfact(int natoms,vector3 atpos[],int G[3])
{
complex S;
int i,j;
real dot;
const real twopi = (real)2.0*M_PI;
S.x = S.y = (real)0.0;
for (i=0; i < natoms; i++) {
dot = (real)0.0;
for (j=0; j < 3; j++)
dot -= atpos[i].v[j]*G[j];
dot *= twopi;
S.x += cos(dot);
S.y += sin(dot);
}
return S;
}
/* Calculates the local pseudopotential in G-space on a parallelopiped
* grid Nx by Ny by Nz in each direction.
* The result is put in Vp in z then y then x order.
* The Nx,Ny,Nz values are in the basis structure, as are the lattice
* vectors.
*
* The pseudopotential in G-space, Vp(G), is given in terms of the
* local pseudopot. in r-space Vp(r) as given below for a single species:
*
* Vp(G) = integral_{unit cell} { d^3r Vp(r)*exp(-i*G*r)/sqrt(Vol) }
* = S(G)*Va(G)/sqrt(Vol)
*
* Here Vp(r) = sum_{R,tau} { Va(r-R-tau) }
* R is a lattice vector and tau is an atomic position in the unit cell.
* Va(r) is the atomic pseudopotential given above and Va(G) is its
* continuous Fourier transform. Vol is the unit cell volume.
*
* ---> Vp(G) HAS sqrt(Vol) INSTEAD OF Vol AS WOULD BE EXPECTED <----
* ---> FOR THE DISCRETE FFT OF Vp(r) <----
*
* A G-vector G_n = n*G where n is a row vector in Z^3 and G is
* G = 2*pi*(latvec)^(-1) where latvec is that 3x3 matrix of primitive
* lattice vectors contained in the columns.
*
* Cutoff issues:
*
* We keep n components within a box satisfying -N_i/2 < n_i < N_i/2.
* The points +/-N_i/2 would cause aliasing problems in the FFT, so we
* keep them out.
*
* G=0 issues:
*
* Va(G=0) is infinite, so we must deal with G=0 in some way.
* Here we choose to set Vp(G=0) = 0...the G=0 energy associated with
* the difference of the pseudopotential and the coulomb potential is
* calculated below.
*
*/
void
Vlocpot(Everything &everything)
{
//24/10/02 ipd: wavelets already have the next three lines
// decide with prof. Arias wheather it should be in basis indep
// place (dft.c? elecvars:setup? ... )
if(everything.ioninfo.read_vloc_flag)
everything.elecvars.Vlocpot.read(everything.ioninfo.vloc_filename);
else{
BasisSpec &basis_spec = everything.basis_spec;
Lattice &lattice = everything.lattice;
Ioninfo &ioninfo = everything.ioninfo;
// Where the result goes!
complex *Vp = everything.elecvars.Vlocpot.data.d;
complex SG;
real G,G2,G2max,s=0.0,svol,r;
real ym1,y0,y1,y2,a,b,c,x;
int beta[3],i,j,index,n[3],Nx,Ny,Nz,sp,nmax=0;
/* Constants */
const real fourpi = (real)4.0*M_PI;
const real sixth = (real)1.0/(real)6.0;
/* Copy Nx, Ny, Nz from basis */
Nx = basis_spec.Nx;
Ny = basis_spec.Ny;
Nz = basis_spec.Nz;
/* Calculate 1/sqrt(unitcellvol) */
svol = (real)1.0/sqrt(lattice.unit_cell_volume);
dft_log("\n----- Vlocpot() -----\n");
dft_log("Calculating local ionic potential with");
dft_log("grid size Nx=%d, Ny=%d, Nz=%d and\n",Nx,Ny,Nz);
dft_log("cell_volume=%lg\n",lattice.unit_cell_volume);
dft_log("latvec = \n");
lattice.latvec.print(dft_global_log,"%lg ");
dft_log("Setting Vpseudo(G=0) = 0\n");
dft_log("\n");
dft_log_flush();
/* Zero out pseudopotential */
for (i=0; i < basis_spec.NxNyNz; i++)
Vp[i].x = Vp[i].y = 0.0;
G2max = (real)0.0;
/* Calculate Vp in G-space: the n sums run over the points we keep in
* the grid (here the FFT G-lattice). */
/* First, loop over the different atomic species */
for (sp = 0; sp < ioninfo.nspecies; sp++) {
Speciesinfo *sp_info = &(ioninfo.species[sp]);
/* Optional report */
dft_log("species = %d Z = %lg natoms = %d\n",
sp, sp_info->Z, sp_info->natoms);
for (i=0; i < sp_info->natoms; i++)
dft_log(DFT_ANAL_LOG,"atpos[%d]=(%10lg,%10lg,%10lg)\n",i,
sp_info->atpos[i].v[0],
sp_info->atpos[i].v[1],
sp_info->atpos[i].v[2]);
dft_log("\n");
dft_log_flush();
nmax = sp_info->pot.ngrid_loc-3;
/* Loop over FFT box and add in contribution from current species */
for (n[0] = -Nx/2+1; n[0] < Nx/2; n[0]++)
for (n[1] = -Ny/2+1; n[1] < Ny/2; n[1]++)
for (n[2] = -Nz/2+1; n[2] < Nz/2; n[2]++) {
if (n[0] < 0) beta[0] = n[0]+Nx;
else beta[0] = n[0];
if (n[1] < 0) beta[1] = n[1]+Ny;
else beta[1] = n[1];
if (n[2] < 0) beta[2] = n[2]+Nz;
else beta[2] = n[2];
index = beta[2]+Nz*(beta[1]+Ny*beta[0]);
G2 = (real)0.0;
for (i=0; i < 3; i++)
for (j=0; j < 3; j++)
G2 += (real)(n[i]*n[j])*lattice.GGT.m[i][j];
if (G2 > G2max)
G2max = G2;
/* G != 0 */
if (n[0] != 0 || n[1] != 0 || n[2] != 0) {
/* Calculate structure factor */
SG = structfact(sp_info->natoms,
sp_info->atpos,n);
/* Interpolate pseudopotential using internal tables */
/* 3rd order fit with j-1, j, j+1, and j+2 points */
G = sqrt(G2);
r = G/sp_info->pot.dq_loc;
j = (int)r;
if (j > nmax)
die("Local pseudopotential grid is too small!|G|=%f\tindex=%d\n", G,j);
if (j < 1)
die("Sampling local pseudopot. too close to 0!\n");
ym1 = sp_info->pot.V_loc[j-1];
y0 = sp_info->pot.V_loc[j ];
y1 = sp_info->pot.V_loc[j+1];
y2 = sp_info->pot.V_loc[j+2];
a = sixth*(y2-ym1+3.0*(y0-y1));
b = 0.5*(y1+ym1-y0-y0);
c = y1-sixth*(y2+3.0*y0+ym1+ym1);
x = r-j;
s = ((a*x+b)*x+c)*x+y0;
/* This line below is quadratic fitting: */
/* s = y0 + 0.5*x*(y1-ym1+x*(y1+ym1-y0-y0)); */
/* Add in -4*pi*Z*e^2/G^2 which is needed */
s -= fourpi*sp_info->Z/G2;
/* Multiply by inverse sqrt(volume) and accumulate */
s *= svol;
Vp[index].x += SG.x*s;
Vp[index].y += SG.y*s;
}
dft_log(DFT_ANAL_LOG,
"n=(%3d,%3d,%3d) G2 = %10.6lf Vp=%15le%+15lei\n",
n[0],n[1],n[2],G2,Vp[index].x,Vp[index].y);
dft_log(DFT_ANAL_LOG,
"SG = %15le%+15lei\n",SG.x,SG.y);
} /* loop on FFT box */
} /* loop on species */
dft_log("Maximum 0.5*|G|^2 = %lg\n\n",0.5*G2max);
dft_log_flush();
}
}
/*
* Fills in dVj with the derivative of Vloc_pseudo versus the j'th coordinate
* of atom 'atom' of species 'species'. Used for calculating forces
* on the atoms.
*/
void
dVlocpot_datom_pos(const Ioninfo &ioninfo,
int species,int atom,
CoeffSpaceScalarFieldColumn &dVx,
CoeffSpaceScalarFieldColumn &dVy,
CoeffSpaceScalarFieldColumn &dVz)
{
BasisSpec &basis_spec = *dVx.basis->basis_spec;
Lattice *lattice = basis_spec.lattice;
real G,G2,vatomic=0.0,svol,r;
real ym1,y0,y1,y2,a,b,c,x;
int beta[3],i,j,index,n[3],Nx,Ny,Nz,nmax=0;
real phase,taux,tauy,tauz;
complex Stau;
real invdq=0.0;
Speciesinfo *sp_info = &(ioninfo.species[species]);
/* Constants */
const real twopi = (real)2.0*M_PI;
const real fourpi = (real)4.0*M_PI;
const real sixth = (real)1.0/(real)6.0;
/* Copy Nx, Ny, Nz from basis */
Nx = basis_spec.Nx;
Ny = basis_spec.Ny;
Nz = basis_spec.Nz;
/* Calculate 1/sqrt(unitcellvol) */
svol = (real)1.0/sqrt(lattice->unit_cell_volume);
/* Zero out the three output arrays */
for (i=0; i < Nx*Ny*Nz; i++)
dVx.data.d[i].x = dVx.data.d[i].y =
dVy.data.d[i].x = dVy.data.d[i].y =
dVz.data.d[i].x = dVz.data.d[i].y = 0.0;
nmax = sp_info->pot.ngrid_loc-3;
invdq = (real)1.0/sp_info->pot.dq_loc;
/* Loop over FFT box */
for (n[0] = -Nx/2+1; n[0] < Nx/2; n[0]++)
for (n[1] = -Ny/2+1; n[1] < Ny/2; n[1]++)
for (n[2] = -Nz/2+1; n[2] < Nz/2; n[2]++) {
if (n[0] < 0) beta[0] = n[0]+Nx;
else beta[0] = n[0];
if (n[1] < 0) beta[1] = n[1]+Ny;
else beta[1] = n[1];
if (n[2] < 0) beta[2] = n[2]+Nz;
else beta[2] = n[2];
index = beta[2]+Nz*(beta[1]+Ny*beta[0]);
G2 = (real)0.0;
for (i=0; i < 3; i++)
for (j=0; j < 3; j++)
G2 += (real)(n[i]*n[j])*lattice->GGT.m[i][j];
/* G != 0 */
if (n[0] != 0 || n[1] != 0 || n[2] != 0) {
/* Interpolate pseudopotential using internal tables */
/* 3rd order fit with j-1, j, j+1, and j+2 points */
G = sqrt(G2);
r = G*invdq;
j = (int)r;
if (j > nmax)
die("Local pseudopotential grid is too small!\n");
if (j < 1)
die("Sampling local pseudopot. too close to 0!\n");
ym1 = sp_info->pot.V_loc[j-1];
y0 = sp_info->pot.V_loc[j ];
y1 = sp_info->pot.V_loc[j+1];
y2 = sp_info->pot.V_loc[j+2];
a = sixth*(y2-ym1+3.0*(y0-y1));
b = 0.5*(y1+ym1-y0-y0);
c = y1-sixth*(y2+3.0*y0+ym1+ym1);
x = r-j;
vatomic = ((a*x+b)*x+c)*x+y0;
/* This line below is quadratic fitting: */
/* s = y0 + 0.5*x*(y1-ym1+x*(y1+ym1-y0-y0)); */
/* Add in -4*pi*Z*e^2/k^2 which is needed */
vatomic -= fourpi*sp_info->Z/G2;
/* Get position of the atom in question and calculate
* Stau = exp(-2*pi*i*dot(G,tau)) */
taux = sp_info->atpos[atom].v[0];
tauy = sp_info->atpos[atom].v[1];
tauz = sp_info->atpos[atom].v[2];
phase = -twopi*(n[0]*taux + n[1]*tauy + n[2]*tauz);
Stau.x = cos(phase);
Stau.y = sin(phase);
/* dV_j = -2*pi*i*n[j]*vatomic/sqrt(Vol)*Stau */
dVx.data.d[index].x = n[0]*twopi*vatomic*svol*Stau.y;
dVx.data.d[index].y = -n[0]*twopi*vatomic*svol*Stau.x;
dVy.data.d[index].x = n[1]*twopi*vatomic*svol*Stau.y;
dVy.data.d[index].y = -n[1]*twopi*vatomic*svol*Stau.x;
dVz.data.d[index].x = n[2]*twopi*vatomic*svol*Stau.y;
dVz.data.d[index].y = -n[2]*twopi*vatomic*svol*Stau.x;
}
} /* loop on FFT box */
}
/*
* Same as above, except Va(|G|), the atomic pseudopotential in
* G-space, is replaced by its derivative versus |G|. The derivatives
* are the exact analytical derivatives of the interpolation algorithm
* in the above code.
*/
void
Vlocpot_prime(Basis *basis,
Lattice *lattice,
Ioninfo *ioninfo,
complex *Vpprime)
{
complex SG;
real G=0.0,G2,G2max,s=0.0,svol,r;
real ym1,y0,y1,y2,a,b,c,x,invdq;
int beta[3],i,j,index,n[3],Nx,Ny,Nz,sp,nmax=0;
/* Constants */
const real eightpi = (real)8.0*M_PI;
const real sixth = (real)1.0/(real)6.0;
/* Copy Nx, Ny, Nz from basis */
Nx = basis->basis_spec->Nx;
Ny = basis->basis_spec->Ny;
Nz = basis->basis_spec->Nz;
/* Calculate 1/sqrt(unitcellvol) */
svol = (real)1.0/sqrt(lattice->unit_cell_volume);
dft_log("\n----- Vlocprime_pseudo() -----\n");
dft_log("Calculating deriv of local pseudopotential with");
dft_log("grid size Nx=%d, Ny=%d, Nz=%d and\n",Nx,Ny,Nz);
dft_log("cell_volume=%lg\n",lattice->unit_cell_volume);
dft_log("latvec = \n");
lattice->latvec.print(dft_global_log,"%lg ");
dft_log("Setting Vpseudo'(G=0) = 0\n");
dft_log("\n");
dft_log_flush();
/* Zero out pseudopotential */
for (i=0; i < basis->basis_spec->NxNyNz; i++)
Vpprime[i].x = Vpprime[i].y = 0.0;
G2max = (real)0.0;
/* Calculate Vp in G-space: the n sums run over the points we keep in
* the grid (here the FFT G-lattice). */
/* First, loop over the different atomic species */
for (sp = 0; sp < ioninfo->nspecies; sp++) {
Speciesinfo *sp_info = &(ioninfo->species[sp]);
/* Optional report */
dft_log("species = %d Z = %lg natoms = %d\n",
sp,sp_info->Z,sp_info->natoms);
for (i=0; i < sp_info->natoms; i++)
dft_log(DFT_ANAL_LOG,
"atpos[%d]=(%10lg,%10lg,%10lg)\n",i,
sp_info->atpos[i].v[0],
sp_info->atpos[i].v[1],
sp_info->atpos[i].v[2]);
dft_log("\n");
dft_log_flush();
nmax = sp_info->pot.ngrid_loc-3;
/* Loop over FFT box and add in contribution from current species */
for (n[0] = -Nx/2+1; n[0] < Nx/2; n[0]++)
for (n[1] = -Ny/2+1; n[1] < Ny/2; n[1]++)
for (n[2] = -Nz/2+1; n[2] < Nz/2; n[2]++) {
if (n[0] < 0) beta[0] = n[0]+Nx;
else beta[0] = n[0];
if (n[1] < 0) beta[1] = n[1]+Ny;
else beta[1] = n[1];
if (n[2] < 0) beta[2] = n[2]+Nz;
else beta[2] = n[2];
index = beta[2]+Nz*(beta[1]+Ny*beta[0]);
G2 = (real)0.0;
for (i=0; i < 3; i++)
for (j=0; j < 3; j++)
G2 += (real)(n[i]*n[j])*lattice->GGT.m[i][j];
if (G2 > G2max)
G2max = G2;
/* G != 0 */
if (n[0] != 0 || n[1] != 0 || n[2] != 0) {
/* Calculate structure factor */
SG = structfact(sp_info->natoms,
sp_info->atpos,n);
/* Interpolate pseudopotential using internal tables */
/* 3rd order fit with j-1, j, j+1, and j+2 points */
G = sqrt(G2);
invdq = (real)1.0/sp_info->pot.dq_loc;
r = G*invdq;
j = (int)r;
if (j > nmax)
die("Local pseudopotential grid is too small!\n");
if (j < 1)
die("Sampling local pseudopot. too close to 0!\n");
ym1 = sp_info->pot.V_loc[j-1];
y0 = sp_info->pot.V_loc[j ];
y1 = sp_info->pot.V_loc[j+1];
y2 = sp_info->pot.V_loc[j+2];
a = sixth*(y2-ym1+3.0*(y0-y1));
b = 0.5*(y1+ym1-y0-y0);
c = y1-sixth*(y2+3.0*y0+ym1+ym1);
x = r-j;
s = ((3.0*a*x+2.0*b)*x+c)*invdq;
/* Add in 8*pi*Z*e^2/k^3 (deriv of coulomb) */
s += eightpi*sp_info->Z/(G2*G);
/* Multiply by inverse sqrt(volume) and accumulate */
s *= svol;
Vpprime[index].x += SG.x*s;
Vpprime[index].y += SG.y*s;
}
dft_log(DFT_ANAL_LOG,
"n=(%3d,%3d,%3d) G2 = %10.6lf Vp'=%15le%+15lei\n",
n[0],n[1],n[2],G2,Vpprime[index].x,Vpprime[index].y);
dft_log(DFT_ANAL_LOG,
"SG = %15le%+15lei\n",SG.x,SG.y);
} /* loop on FFT box */
} /* loop on species */
dft_log("Maximum 0.5*|G|^2 = %lg\n\n",0.5*G2max);
}
/*
* The G=0 contribution to the energy for the local pseudopotentials.
* For each species, this is equal to
*
* integral_{unit cell} { d^3r n(r)*(Vp(r)-Vc(r))|G=0 }
*
* where n(r) is the number density of the electrons, Vp(r) is the
* pseudopotential (periodic), and Vc(r) is the potential due to
* a constant density of positive charge equal to the total ionic charges;
* i.e., natoms*Z. The |G=0 denotes that we only need the G=0 part of the
* discrete fourier expansion of Vp(r)-Vc(r).
*
* Looking at the above routine's comments, the discrete fourier transform
* of Vp(r) is Vp(G) = S(G)*Va(G)/Vol, where Vol=unit cell volume.
* Similarly, Vc(G) = -S(G)*4*pi*Z*e^2/(G^2*Vol).
* Since S(0) = natoms = nions,
*
* (Vp(r)-Vc(r))|G=0 = (natoms*lim_{G->0} {Va(G)+4*pi*Z*e^2/G^2})/Vol,
* which is a constant. The integral of n(r) then gives nelectrons
* and the final result is nelectrons*(Vp(r)-Vc(r))|G=0.
*
* report <= 0: silent
* report > 0: says what it is doing
*
*/
real
Vlocpot_GzeroEnergy(real nelectrons, const Lattice &lattice,
const Ioninfo &ioninfo)
{
real e,etot;
int sp;
dft_log("\n----- VpseudoGzeroEnergy() -----\n");
dft_log("Nelectrons = %lg unit cell volume = %lg\n",
nelectrons,lattice.unit_cell_volume);
e = etot = 0.0;
for (sp = 0; sp < ioninfo.nspecies; sp++) {
Speciesinfo *sp_info = &(ioninfo.species[sp]);
e = (real)(nelectrons* sp_info->natoms)*
sp_info->pot.V_loc[0]/lattice.unit_cell_volume;
dft_log("species = %d Vloc[0] = %19.12le e=%19.12le\n",
sp, sp_info->pot.V_loc[0], e);
etot += e;
}
dft_log("\nCore energy = %19.12le\n",etot);
return etot;
}
/*
* the routine below calculates (and prints) the charge density
* at the nucleus of the ions.
*
* rho_k(r) = sum_n f(n) |sum_G C_nk(G) exp(i(G*r))|^2
*
*/
void calc_n_ion(Everything &e)
{
// column_bundle *Y = e.elecvars.Y;
Ioninfo &ioninfo = e.ioninfo;
Elecinfo &einfo = e.elecinfo;
// Basis *basis = e.basis;
#ifdef SCALAR_IS_COMPLEX
vector3 G, r;
complex i(0,1);
real scl = 1.0/e.lattice.unit_cell_volume;
dft_log("\n >>>>>>>>> Charge density at the ion positions <<<<<<<<<<\n\n");
for(int sp = 0; sp < ioninfo.nspecies; sp++)
for(int ion = 0; ion < ioninfo.species[sp].natoms; ion++){
r = e.lattice.R*ioninfo.species[sp].atpos[ion];
for(int q = 0; q < einfo.nstates; q++){
BlochState &s = e.elecvars.states[q];
real rho = 0;
for(int b = 0; b < einfo.nbands; b++){
complex tmp(0,0);
for(int n = 0; n < s.basis.nbasis; n++){
G.v[0] = s.basis.Gx[n];
G.v[1] = s.basis.Gy[n];
G.v[2] = s.basis.Gz[n];
G = e.lattice.GT*G;
tmp += s.Y.col[b]->data.d[n]*exp(i*(G*r));
}
rho += scl*s.w*REAL(s.F.c[b])*abs2(tmp);
}
dft_log("q = %3d, rho(%d.%d:%25.16e) = %25.16e\n", q, sp, ion, G*r, rho);
dft_log_flush();
}
}
dft_log("\n>>>>>>>>>>>>>>> done <<<<<<<<<<<<<<<<<<\n\n");
#elif defined SCALAR_IS_REAL
#error "Don't know how to do manual fft for real in calc_n_ion()"
#else
#error scalar is neither complex nor real
#endif
}
syntax highlighted by Code2HTML, v. 0.9.1