/*
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 the U matrices, the C column_bundles, and charge density n.
*
*/
/* $Id: calcUVCn.cpp,v 1.13.2.12 2003/05/29 18:54:13 ivan Exp $ */
#include "header.h"
/*
* Calculate the U matrices and their various powers:
*
* U[k] = Y[k]^O(Y[k])
*/
void calc_U(Elecinfo &einfo,Elecvars &evars)
{
#ifdef DFT_PROFILING
timerOn(14); // Turn on calc_U timer
#endif // DFT_PROFILING
int q;
for (q=0; q < einfo.nstates; q++) {
BlochState &s=evars.states[q];
ColumnBundle OY(s.Y);
/* does: U[k] = Y[k]^(O(Y[k])); */
O(s.Y,OY);
s.U = s.Y^OY;
s.U.hermetian = 1; // U is hermetian, so set the flag
s.Umhalf = Uminusonehalf(s.U,s.W,s.mu);
}
#ifdef DFT_PROFILING
timerOff(14); // Turn off calc_U timer
#endif // DFT_PROFILING
}
/*
* Calculate the subspace rotation matrix : V[q] = exp(iB[q])
*/
void
calc_V(Elecinfo &einfo, Elecvars &evars)
{
int q,i;
/* If we aren't doing subspace rotations, then we should do nothing! */
if (!einfo.subspace_rotation)
return;
#ifdef DFT_PROFILING
timerOn(41); // Calc V timer
#endif
diag_matrix expibeta(evars.states[0].B.nr);//assuming everybody is the same
// consider moving it inside the loop
for (q=0; q < einfo.nstates; q++) {
BlochState &s=evars.states[q];
/* B[k] = Z[k]*beta[k]*Zdag[k] */
diagonalize_herm(s.beta,s.Z,s.B,s.B.nr);
/* V[k] = Z[k]*exp(i*beta[k])*Zdag[k] */
for (i=0; i < s.B.nr; i++) {
#ifdef SCALAR_IS_COMPLEX
complex ibeta;
ibeta.x = 0.0;
ibeta.y = s.beta[i];
expibeta.c[i] = exp(ibeta);
#elif defined SCALAR_IS_REAL
#error "Don't know how to do subspace rotations for real"
#else
#error scalar is neither complex nor real
#endif
}
s.V = s.Z*(expibeta*herm_adjoint(s.Z));
}
#ifdef DFT_PROFILING
timerOff(41); // Calc V timer
#endif
}
/*
* Calculate the orthonormal wavefunctions:
*
* state[q]: C = Y*Umhalf*Vdag
*/
void
calc_C(Elecinfo &einfo,Elecvars &evars)
{
#ifdef DFT_PROFILING
timerOn(15); // Turn on calc_C timer
#endif // DFT_PROFILING
int q;
//just loop over the bloch states and transform Y to C
for (q=0; q < einfo.nstates; q++) {
BlochState &s = evars.states[q];
if (einfo.subspace_rotation)
/* C = Y*Umhalf; */
do_ColumnBundle_matrix_mult(s.Y, s.Umhalf*herm_adjoint(s.V), s.C, 0);
else
do_ColumnBundle_matrix_mult(s.Y, s.Umhalf,s.C, 0);
}
#ifdef DFT_PROFILING
timerOff(15); // Turn off calc_C timer
#endif // DFT_PROFILING
}
/* Combines U, V, and C */
void
calc_UVC(Elecinfo &einfo,Elecvars &evars)
{
calc_U(einfo,evars);
calc_V(einfo,evars);
calc_C(einfo,evars);
}
/* Calculate electron density: n = sum_k { w[q]*diagouter(F[q],I*C[q]) } */
void
calc_n(Elecinfo &einfo, Elecvars &evars,
Lattice &lattice, Symmetries &symm)
{
#ifdef DFT_PROFILING
timerOn(16); // Turn on calc_n timer
#endif // DFT_PROFILING
/* Get local pointer to the relevant variables to make the formulae
* legible */
RealSpaceScalarFieldColumn &n = evars.n; /* reference */
RealSpaceScalarFieldColumn &n_up = evars.n_up; /* reference */
RealSpaceScalarFieldColumn &n_dn = evars.n_dn; /* reference */
// RealSpaceScalarFieldColumn &n_ud = evars.n_ud; /* reference */
// RealSpaceScalarFieldColumn &n_du = evars.n_du; /* reference */
//ipd: these guys are always present, no need of if sttatement
n.zero_out();
n_up.zero_out();
n_dn.zero_out();
//ipd: change the += below w/ multadds to save space
for (int q=0; q < einfo.nstates; q++){
BlochState &s = evars.states[q];
switch(einfo.spintype){
case NOSPIN: // NOSPIN MODE
n += ((scalar)s.w)*diagouterI(s.F,s.C); break;
case ZSPIN: // ZSPIN MODE
switch (s.qnum.spin){
case 1 : n_up += ((scalar)s.w)*diagouterI(s.F,s.C); break;
case -1 : n_dn += ((scalar)s.w)*diagouterI(s.F,s.C); break;
default : die("Bad spin designation in ZSPIN mode: %d", s.qnum.spin);
}
break;
default:
die("Invalid spintype: %d", einfo.spintype);
}
}
switch(einfo.spintype){
case NOSPIN:
// symmetrize the total charge density.
symmetrize_n(n, &symm);
break;
case ZSPIN:
// symmetrize
symmetrize_n(n_up, &symm);
symmetrize_n(n_dn, &symm);
n = n_up; n+= n_dn;
break;
default:
die("Invalid spintype: %d", einfo.spintype);
}
// integral of n should be number of electrons per unit cell.
/*ipd: it's more complicated for wavelets
-- several choices:
* do s^(JN) - may be expensive
* add just the top level grids - terrible
* make a basis dependent function still have to deal w/ it's cost
real sum = REAL(sum_vector(n));
sum *= lattice.unit_cell_volume / n.n;
if ( fabs(sum/einfo.nelectrons-1) > 1e-4 )
{
dft_log(DFT_ANAL_LOG,
">>>>Real space charge integrates to %g != %g (nelectrons)\n",
sum, einfo.nelectrons);
die("Quiting due to violation of charge conservation.");
}
*/
#ifdef DFT_PROFILING
timerOff(16); // Turn off calc_n timer
#endif // DFT_PROFILING
}
// Combines U, V, C and n
void calc_UVCn(Elecinfo &einfo,Elecvars &evars,
Lattice &lattice,Symmetries &symm)
{
calc_U(einfo,evars);
calc_V(einfo,evars);
calc_C(einfo,evars);
calc_n(einfo,evars,lattice,symm);
}
syntax highlighted by Code2HTML, v. 0.9.1