/*
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.
*/
/*
* Sohrab Ismail-Beigi June 12, 1997
*
* A set of routines to calculate Fermi-Dirac fillings for the
* electronic bands.
*
*/
/* $Id: fermifill.cpp,v 1.10.2.3 2003/05/29 18:54:24 ivan Exp $ */
#include "header.h"
// #include "parallel.h"
/* The Fermi-Dirac distribution */
static real
Fermi_Dirac_filling(const SpinType spintype,const real epsilon,
const real kT,const real mu)
{
real ex = (epsilon-mu)/kT;
if (ex < 700) { // prevent floating point number overflow
if (spintype == NOSPIN)
return 2.0/(1.0+exp(ex));
else if (spintype == ZSPIN)
return 1.0/(1.0+exp(ex));
else
die("Spintype not supported in Fermi_Dirac_filling.\n");
}
return 0.0;
}
/*
* Calculates the number of electrons given kT and mu and the
* the eigenenergies epsilon[k][band] using the Fermi-Dirac distribution
*/
static real
Nelecs_fermi(const SpinType spintype,
const int nstates,const int nbands,
real*w, real **epsilon,const real kT,const real mu)
{
int q;
int band;
real N;
N = (real)0.0;
for (q=0; q < nstates; q++)
for (band=0; band < nbands; band++)
N += w[q]*Fermi_Dirac_filling(spintype,epsilon[q][band],kT,mu);
return N;
}
/*
* The function calculates the fillings of the bands at all kpts by adjusting
* the chemical potential mu (via binary sections) so as to make
* N(mu) = elecinfo.nelectrons. Then it calculates the fillings based on
* that mu.
*
* kT is the temperature in Hartrees.
*
* The eigenenergies must be sorted in ascending order and be placed
* in elecvars.Hsub_eigs[k][band].
*
*/
void
calc_fermi_fillings(Everything &everything)
{
Elecinfo &elecinfo = everything.elecinfo;
Elecvars &elecvars = everything.elecvars;
Control &cntrl = everything.cntrl;
real kT = elecinfo.kT;
int nstates = elecinfo.nstates;
int nbands = elecinfo.nbands;
int q,band;
real mu1,mu2,mu;
real N;
const real Nideal = elecinfo.nelectrons;
real **epsilon; // contains the eigenvalues from each BlochState
real *w; // contains the weights from each BlochState
int section;
// These two arrays are needed for Nelecs_fermi, and provide a convenient
// shorthand.
// Collect eigenvalues to create the epsilon array.
epsilon = (real **) mymalloc(sizeof(real*)*nstates,"epsilon",
"calc_fermi_fillings()");
for(q=0; q<nstates; q++)
epsilon[q] = elecvars.states[q].Hsub_eigs;
// epsilon[q][b] gives the eigenvalue for band b in state q
// Collect weights to create the w array.
w = (real*) mymalloc(sizeof(real)*nstates,"w",
"calc_fermi_fillings()");
for(q=0; q<nstates; q++)
w[q] = elecvars.states[q].w;
// w[q] gives the weight for state q
#define NSECTIONS 50
/* Set mu1 = minimal eigenenergy, mu2=maximal eigenenergy */
mu1 = mu2 = epsilon[0][0];
for (q=0; q < nstates; q++)
for (band=0; band < nbands; band++)
{
if (epsilon[q][band] < mu1) mu1 = epsilon[q][band];
if (epsilon[q][band] > mu2) mu2 = epsilon[q][band];
}
dft_log("\n--- calc_fermi_fillings() ---\n");
dft_log("nstates = %d nbands=%d kT = %le Nideal=%le\n",
nstates,nbands,kT,Nideal);
dft_log("min(epsilon)=%le max(epsilon)=%le\n",mu1,mu2);
dft_log("states and eigenenergies follow:\n\n");
for (q=0; q < nstates; q++)
{
dft_log("q_k[%d] = [ %lg %lg %lg ]\n",q,
elecvars.states[q].qnum.kvec.v[0],
elecvars.states[q].qnum.kvec.v[1],
elecvars.states[q].qnum.kvec.v[2]);
for (band=0; band < nbands; band++)
dft_log("%le\t",epsilon[q][band]);
dft_log("\n\n");
}
/* Do binary sections NSECTIONS times to locate mu so N(mu) = Nideal */
for (section = 0; section < NSECTIONS; section++)
{
mu = (mu1 + mu2)/2.0;
N = Nelecs_fermi(elecinfo.spintype,nstates,nbands,w,epsilon,kT,mu);
dft_log("section=%3d mu=[%le,%le,%le] N=%le\n",
section,mu1,mu,mu2,N);
// if converged, break out
if (fabs(N - Nideal) < cntrl.nel_fermi_tol) break;
// bisection search:
if ( N > Nideal )
mu2 = mu;
else
mu1 = mu;
}
// If we couldn't converge, die (this should never happen
// as 2^NSECTIONS should reprsent a ridiculously small tolerance...
if (section == NSECTIONS)
die("\nCould not find appropriate mu in %d bisections!\nQuitting\n",
NSECTIONS);
/* Store mu into elecinfo */
elecinfo.mu = mu;
/* Print final mu and N */
dft_log("\nFinal mu=%le N=%le N-Nideal=%le\n\n",mu,N,N-Nideal);
/* Calculate the fillings for the mu we've found and print out fillings */
dft_log("Calculated states, weights, and fillings follow:\n\n");
for (q=0; q < nstates; q++)
{
dft_log("%15.10lf %15.10lf %15.10lf %15.9le\n",
elecvars.states[q].qnum.kvec.v[0],
elecvars.states[q].qnum.kvec.v[1],
elecvars.states[q].qnum.kvec.v[2],
elecvars.states[q].w);
for (band=0; band < nbands; band++)
{
elecvars.states[q].F.c[band] =
Fermi_Dirac_filling(elecinfo.spintype,epsilon[q][band],kT,mu);
dft_log("%15.9le ",REAL(elecvars.states[q].F.c[band]));
}
dft_log("\n#\n");
}
dft_log("\n");
// Free the epsilon array which we created at the beginning.
// Note: This doesn't affect the Hsub_eigs pointers contained in epsilon,
// does it???
myfree(epsilon);
myfree(w);
}
syntax highlighted by Code2HTML, v. 0.9.1