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