/*
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 Mar. 29, 1997
*
* Calculate all the various energy terms
*
*/
/* $Id: calcener.cpp,v 1.11.2.26 2003/05/29 18:54:14 ivan Exp $ */
#include "header.h"
/*
* Calclate kinetic energy:
*
* KE = -0.5*sum_k { w[k]*trace(diag(C[k]^F[k]*L(C[k]))) }
*
*/
void calc_KE(Everything &everything)
{
#ifdef DFT_PROFILING
timerOn(18); // Turn on calc_KE timer
#endif // DFT_PROFILING
Elecinfo &einfo = everything.elecinfo;
Elecvars &evars = everything.elecvars;
Energies &ener = everything.energies;
int q; // quantum state loop counter
/* Create temporary workspace LC*/
// now loop over states to accumulate the kinetic energy
ener.KE = 0.0;
for (q=0; q < einfo.nstates; q++){
BlochState &s = evars.states[q];
//ipd: for now create a copy - later use scratch space
ColumnBundle LC(s.C);
// ipd: ColumnBundle is much simpler than column_bundle
// therefore we do not need to coppy the innards
// copy_innards_column_bundle(s.C,&LC);
/* does: ener.KE += REAL(sum_vector(diaginner(F[q],C[q],L(C[q]))))*w[q]; */
L(s.C,LC);
ener.KE += REAL(trace(diaginner(s.F,s.C,LC)))*s.w;
}
ener.KE *= -0.5;
#ifdef DFT_PROFILING
timerOff(18); // Turn off calc_KE timer
#endif // DFT_PROFILING
}
/*
* Hartree energy.
*
* The routine assumes that the poisson equation has already been solved
* and uses the electrostatic potential in evars.d.
*
*/
void calc_EH(Everything &everything)
{
#ifdef DFT_PROFILING
timerOn(21); // Turn on calc_EH timer
#endif // DFT_PROFILING
Elecvars &evars = everything.elecvars;
RealSpaceScalarFieldColumn &n = evars.n;
CoeffSpaceScalarFieldColumn &d = evars.d;
//ipd: for WL the "^" operator changes n, just because
//of our realspace conventions where we keep the ghosts
//matched - try not to use it
//everything.energies.EH = 0.5*REAL(n^Jdag(Obar(d)));
everything.energies.EH = 0.5*REAL(J(n)^(Obar(d)));
#ifdef DFT_PROFILING
timerOff(21); // Turn off calc_EH timer
#endif // DFT_PROFILING
}
// Exchange correlation energy
void calc_Exc(Everything &everything)
{
#ifdef DFT_PROFILING
timerOn(22); // Turn on calc_Exc timer
#endif // DFT_PROFILING
Elecinfo &einfo = everything.elecinfo;
Elecvars &evars = everything.elecvars;
Energies &ener = everything.energies;
if (einfo.spintype == FREESPIN)
die("FREESPIN is not implemented.");
RealSpaceScalarFieldColumn &n = evars.n;
switch (einfo.exc_opt) {
case DFT_EXCORR_GGA_PBE: {
/* Evaluate gradient on real space grid */
RealSpaceScalarFieldColumn Dx(n),Dy(n),Dz(n);
int dir,dag;
apply_D(dir=0,dag=0,n,Dx);
apply_D(dir=1,dag=0,n,Dy);
apply_D(dir=2,dag=0,n,Dz);
ener.Exc = REAL(J(n)^(O(J(exc_GGA(n,Dx,Dy,Dz)))));
}
break;
case DFT_EXCORR_LDA: {
/* save on temps:
exc(n,foo);
J(foo, foo);
O(foo,foo);
Jdag(foo,foo);
*/
ener.Exc = REAL(J(n)^(O(J(exc_LDA(n)))));
}
break;
// For now, PW exchange-correlation is our only option. --mhe
/*
case DFT_EXCORR_GGA_PW91:
if (einfo.spintype == NOSPIN)
ener.Exc = REAL(n^Jdag(O(J(exc_LDA(n)+ex_GC(n)))));
else
die("GGA_PW91 not implemented with spin\n");
break;
*/
case DFT_EXCORR_LSD_TETER:
{
if (einfo.spintype == NOSPIN)
die("Can't use LSD_TETER without spin\n");
RealSpaceScalarFieldColumn &n_up = evars.n_up;
RealSpaceScalarFieldColumn &n_dn = evars.n_dn;
RealSpaceScalarFieldColumn exc(n);
exc = exc_LSD_TETER(n_up, n_dn);
ener.Exc = REAL(J(n)^(O(J(exc))));
break;
}
default:
{
die("This exchange correlation (%d) is not supported.",einfo.exc_opt);
break;
}
}
#ifdef DFT_PROFILING
timerOff(22); // Turn off calc_Exc timer
#endif // DFT_PROFILING
}
// Ewald energy
void calc_Eewald(Everything &everything)
{
#ifdef DFT_PROFILING
timerOn(24); // Turn on calc_Eewald timer
#endif // DFT_PROFILING
everything.energies.Eewald = Ewald(everything.ioninfo,
everything.lattice);
#ifdef DFT_PROFILING
timerOff(24); // Turn on calc_Eewald timer
#endif // DFT_PROFILING
}
/*
* Calculate local pseudopotential energy:
*
* Eloc = real{ n^Jdag(Vlocps) }
*
* THE ROUTINE USES THE LOCAL PSEUDOPOTENTIAL IN evars.Vlocps.
* It does NOT recalculate evars.Vlocps.
*
*/
void calc_Eloc(Everything &everything)
{
#ifdef DFT_PROFILING
timerOn(19); // Turn on calc_Eloc timer
#endif // DFT_PROFILING
Elecvars &evars = everything.elecvars;
RealSpaceScalarFieldColumn &n = evars.n;
CoeffSpaceScalarFieldColumn &Vlocpot = evars.Vlocpot;
everything.energies.Eloc = REAL(J(n)^(Vlocpot) );
#ifdef DFT_PROFILING
timerOff(19); // Turn off calc_Eloc timer
#endif // DFT_PROFILING
}
/* Total sum of all the above: sets ener.Etot to the sum. */
void calc_Etot(Energies &ener)
{
/* TAA ener.Etot = ener.Exc; */
ener.Etot =
ener.KE +
ener.Eloc +
ener.Enl +
ener.EH +
ener.Exc +
ener.Ecore +
ener.Eewald +
ener.Epulay;
}
/*
* Calculates the Helmholtz free energy for the system : F = Etot - T*S
* Etot is assumed to be already calculated (ener.Etot). The entropy S,
* the temperature kT, the chemical potential mu, and the fillings and weights
* are all inside of einfo. This is used only in the case were we calculate
* fermi fillings (i.e. einfo.calc__fillings == 1).
*
* If w[q] are the state weights and f[q][i] is the i'th band's filling
* at the q'th state, then (here fillings obey 0<f[q][i]<1 and
* we assume we have a spin-compensated system, spintype = NOSPIN)
*
* T*S = -2.0*sum_q { w[q]*sum_i { f[q][i]*ln(f[q][i]) +
* (1-f[q][i])*ln(1-f[q][i]) } }
*
*/
void calc_F(Everything &e)
{
real E,TS,fqi,kT;
int q,i;
if (e.elecinfo.calc_fillings == 0){
// Don't die, just set F = Etot (==> S = 0, for constant fillings).
// die("\n\ncalc_F() called with einfo.calcfillings == 0!!!\n\n");
dft_log("For constant fillings, F = Etot.");
e.energies.F = e.energies.Etot;
return;
}
// Calculate T*S, and N
TS = (real)0.0;
kT = e.elecinfo.kT;
E = e.energies.Etot;
for (q=0; q < e.elecinfo.nstates; q++)
for (i=0; i < e.elecinfo.nbands; i++){
fqi = REAL(e.elecvars.states[q].F.c[i])/2.0;
if (fqi < 0.0 || fqi > 1.0)
die("\n\nf[%d][%d] is not in [0.0,1.0] in calc_G()!\n\n", q,i);
// If fqi is not 0.0 or 1.0 (to one part in 10^30), then add its
// contribution to TS
if ( fqi > 1.0e-30 && (1.0-fqi) > 1.0e-30)
TS -= kT*2.0*e.elecvars.states[q].w*
( fqi*log(fqi) + (1.0-fqi)*log(1.0-fqi) );
}
e.energies.F = E - TS;
}
/*
* Calculate the energy terms depenent on the electronic variables,
* and recalculate the total energy.
*/
void calc_elec_dependent_energies(Everything &everything)
{
Energies &ener = everything.energies;
Elecinfo &einfo = everything.elecinfo;
calc_KE(everything);
calc_Eloc(everything);
calc_Enl(everything);
calc_EH(everything);
calc_Exc(everything);
calc_Etot(ener); //add all the pieces
if (einfo.calc_fillings == 1)
calc_F(everything);
}
/*
* Calculate the core and Ewald energies only and recalculate the total.
*/
//calc_ionic_energies(Everything &everything); //basis dependent
void calc_core_ewald_pulay_energies(Everything &everything)
{
Energies &ener = everything.energies;
Elecinfo &elecinfo = everything.elecinfo;
calc_Ecore(everything);
calc_Eewald(everything);
calc_Epulay(everything);
calc_Etot(ener);
if (elecinfo.calc_fillings == 1)
calc_F(everything);
}
/*
* Calculates ALL of the above energies
*/
//basis indep
void calc_all_energies(Everything &everything)
{
calc_elec_dependent_energies(everything);
calc_core_ewald_pulay_energies(everything);
}
syntax highlighted by Code2HTML, v. 0.9.1