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