/*
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.
*/
#include "header.h"
#include <time.h>
//
// The following four routines are just to make the code
// further below easier to understand: they just do nothing
// if subrot is off (subspace rotation), but do the requisite
// activity if on.
//
// Calculate square length only if subrot is on
static real
abs2(int nmat,Matrix **M,int subrot)
{
if (subrot)
return abs2(nmat,M);
else
return 0.0;
}
// Calculate dot product only if subrot is on
static real
dot(int nmat,Matrix **M1,Matrix **M2,int subrot)
{
if (subrot)
return REAL(dot(nmat,M1,M2));
else
return 0.0;
}
// do the sum only if subrot is on
static void
scaled_sum(int n,scalar s1,Matrix **m1,scalar s2,Matrix **m2,
Matrix **mout,int subrot)
{
if (subrot)
scaled_sum(n,s1,m1,s2,m2,mout);
}
// do the scaled accumulation only if subrot is on
static void
scale_accumulate(int n,scalar s,Matrix **min,Matrix **mout,int subrot)
{
if (subrot)
scale_accumulate(n,s,min,mout);
}
//
// Calculates U, V, C, n, then solves poisson's equation (d).
// This is the bare minimum of what we must do to get a working
// electronic state.
//
void
calc_UVCn_d(Everything &everything)
{
dft_log_flush();
calc_UVCn(everything.elecinfo,
everything.elecvars,
everything.lattice,
everything.symm);
solve_poisson(everything.elecvars.n,everything.elecvars.d);
}
//
// Calculates U, V, C, n, then solves poisson's equation (d), and
// then calculates the energies that depend on the electronic variables
// (i.e. not the "core" energy or the Ewald sums).
//
void
calc_UVCn_d_elec_dependent_energies(Everything &everything)
{
#ifdef DFT_PROFILING
timerOn(11); // Turn on calc_UCn_d_elec_dependent_e timer
#endif // DFT_PROFILING
calc_UVCn_d(everything);
dft_log("Computing electronic energies\n");
dft_log_flush();
calc_elec_dependent_energies(everything);
#ifdef DFT_PROFILING
timerOff(11); // Turn off calc_UCn_d_elec_dependent_e timer
#endif // DFT_PROFILING
}
//
// Calculates U,V,C,n,d,Hsub,and fillings: this is useful when
// we want to recompute Fermi fillings, where we don't need the
// energies, but do need the electronics state and the subspace
// Hamiltonian and its eigenvalues.
//
void
calc_UVCn_d_Hsub_fillings(Everything &everything)
{
calc_UVCn_d(everything);
dft_log("Computing Hsub\n");
dft_log_flush();
calc_Hsub(everything);
calc_fermi_fillings(everything);
}
//
// Calculate current energy and possibly recalculate fillings if
// the iteration number is a modulus of the nubmer to recompute.
// If iter==-1, then the fillings are not recomputed.
//
static real
calc_electronic_energy(Everything &e,int iter)
{
// If it is time, recalculate the Fermi fillings
if (e.elecinfo.calc_fillings == 1)
if (iter!=-1 && iter%e.elecinfo.niter_recalc_fillings==0)
calc_UVCn_d_Hsub_fillings(e);
// Calculate the energy
calc_UVCn_d_elec_dependent_energies(e);
// Return energy
return e.energies.Etot;
}
//
// Calculate current band struture energy:
//
// E_band = sum_q { tr(C[q]^H*C[q]) }
//
static real
calc_band_struct_energy(Everything &e,int iter)
{
// Orthonormalize wave functions, but *DO NOT* recompute
// electron density (n)
calc_UVC(e.elecinfo,e.elecvars);
if(e.elecinfo.read_n_flag == TRUE) // read the density
solve_poisson(e.elecvars.n, e.elecvars.d);
calc_Hsub(e);
e.energies.Etot = 0;
for (int q=0; q < e.elecinfo.nstates; q++)
e.energies.Etot += REAL( trace(e.elecvars.states[q].Hsub) );
// Return energy
return e.energies.Etot;
}
//
// Print a banner with iteration number, time and date, energies,
// and the band structure.
//
static void
print_electronic_energy(Everything &e, int iter)
{
// Time stamp and print out current energies and subspace Hamiltonian
time_t timenow = time(0);
dft_log("------------------------------------------------------\n");
dft_log("Iteration %d %s\n",iter,ctime(&timenow));
e.energies.print(dft_global_log);
dft_log("\n");
print_Hsub_eigs(e,dft_global_log);
dft_log("\n");
dft_log_flush();
}
//
// Print a banner with iteration number, time and date, band
// energy and the energy eigenvalues.
//
static void
print_band_struct_energy(Everything &e, int iter)
{
// Time stamp and print out current energies and subspace Hamiltonian
time_t timenow = time(0);
dft_log("------------------------------------------------------\n");
dft_log("Iteration %d %s\n",iter,ctime(&timenow));
dft_log("E_band = %21.13le\n\n",e.energies.Etot);
print_Hsub_eigs(e,dft_global_log);
dft_log("\n");
dft_log_flush();
}
//
// Does a line minimization along the search diretion (Ydir/Bdir),
// and returns what it finds to be the
// optimum move along the search direction.
//
// The alogithm used does a quadratic fit by using the current
// energy and directional derivative along the serach dir. plus
// the value of the energy when moved by a stepsize along the search dir.
// The routine then moves to the minimum of the quadratic fit.
//
// The routine has some "error catching" in that it will not accept
// a proposed move if the motion along the direction is either too small
// or too large compared to the serach step size; in such a case, it
// either decreases or increases the step size and tries again. It
// will try only a finite number of times (for now fixed by the #define
// below) before giving up. Also, there is a minimum size of a step
// below which it will refuse to go (to avoid paranoid situations or
// when the noise in the energy function makes minimization hard).
//
static real
linmin(ColumnBundle **Ydir,ColumnBundle **Ygrad,
int subrot,Matrix **Bdir, Matrix **Bgrad,
real (*calc_energy)(Everything &e,int iter),
Everything &e)
{
#ifdef DFT_PROFILING
timerOn(9); // Turn on do_linmin timer
#endif // DFT_PROFILING
// Get the current electronic variables
int nstates = e.elecinfo.nstates;
ColumnBundle **Y = e.elecvars.Y;
Matrix **B = e.elecvars.B;
// Some useful values
real stepsize = e.cntrl.elec_stepsize;
real gamma,dderiv,curvature;
// Save initial energy
const Energies ener0 = e.energies;
// Directional derivative along dir
dderiv = 2.0*dot(nstates,Ygrad,Ydir) + dot(nstates,Bgrad,Bdir,subrot);
// Try to do line minimzations until you succeed or the limit is reached
#define DO_LINMIN_LIMIT 5
for (int i = 0; i < DO_LINMIN_LIMIT; i++) {
// Shift Y by stepsize: Y[q] += stepsize*Ydir[q]
scale_accumulate(nstates,stepsize,Ydir,Y);
//IPD: no scale_accumulate(1,2,3,4,**5**)???
scale_accumulate(nstates,stepsize,Bdir,B,subrot);
// Get the energy
calc_energy(e,-1);
// Shift back
scale_accumulate(nstates,-stepsize,Ydir,Y);
scale_accumulate(nstates,-stepsize,Bdir,B,subrot);
// Do a parabolic fit to the Etot0, Etot, and dderiv
// to get curvature for quadratic fit and proposed minimum (gamma)
curvature = 2.0*(e.energies.Etot-ener0.Etot-stepsize*dderiv) /
(stepsize*stepsize);
gamma = -dderiv/curvature;
dft_log("E0 = %25.15le E = %25.15le\n",ener0.Etot,e.energies.Etot);
dft_log("dderiv = %14.7le curvature = %14.7le\n",dderiv,curvature);
dft_log("stepsize = %14.7le gamma = %14.7le\n",stepsize,gamma);
dft_log_flush();
// If curvature is wrong way, take a bigger step
if (curvature < 0.0) {
stepsize *= 4.0;
continue;
}
// If the proposed minimum is much larger than the stepsize,
// increase stepsize
else if (fabs(gamma/stepsize) > 10.0) {
stepsize *= 10.0;
continue;
}
// If much smaller, decrease
// Tairan (10/15/1998): but the stepsize should not be too small.
// otherwise, numerical error will overwhelm the calculation
// Y += stepsize * dir
else if ( fabs(gamma/stepsize)<0.1 &&
fabs(stepsize)>e.cntrl.stepsize_min ) {
stepsize *= 0.1;
continue;
}
// Otherwise, it was a good linmin so stop!
else
break;
}
#ifdef DFT_PROFILING
timerOff(9); // Turn off do_linmin timer
#endif // DFT_PROFILING
// If update_stepsize flag is set, use gamma as the new stepsize
// if it is larger than the stepsize minimum.
// Otherwise, update stepsize only if it has been updated by
// factors of 10 or 0.1
if (e.cntrl.update_stepsize && fabs(gamma)>e.cntrl.stepsize_min) {
dft_log("stepsize updated to %e : (min is %e)\n",
gamma, e.cntrl.stepsize_min);
dft_log_flush();
e.cntrl.elec_stepsize = gamma;
}
// Return optimum step found
return gamma;
}
//
// Does a line minimization along the search diretion (Ydir/Bdir),
// and returns what it finds to be the
// optimum move along the search direction.
//
// This version uses only derivatives, instead of the energy.
// this is slower, but we can converge the wavefunctions better.
//
// it makes a default step, computes the new directional derivative,
// then extrapolates linearly to where the dderiv would be zero.
//
// The routine has some "error catching" in that it will not accept
// a proposed move if the motion along the direction is either too small
// or too large compared to the serach step size; in such a case, it
// either decreases or increases the step size and tries again. It
// will try only a finite number of times (for now fixed by the #define
// below) before giving up. Also, there is a minimum size of a step
// below which it will refuse to go (to avoid paranoid situations or
// when the noise in the energy function makes minimization hard).
//
static real
linmin_derivs_only(ColumnBundle **Ydir, ColumnBundle **Ygrad,
int subrot, Matrix **Bdir, Matrix **Bgrad,
void (*calc_grad_pgrad)(ColumnBundle **Ygrad, ColumnBundle **pYgrad,
Matrix **Bgrad, Matrix **pBgrad,
int precond, Everything &e),
Everything &e)
{
dft_log("Doing linmin with derivatives only\n");
dft_log_flush();
#ifdef DFT_PROFILING
timerOn(9); // Turn on do_linmin timer
#endif // DFT_PROFILING
// Get the current electronic variables
int nstates = e.elecinfo.nstates;
ColumnBundle **Y = e.elecvars.Y;
Matrix **B = e.elecvars.B;
// Some useful values
real stepsize = e.cntrl.elec_stepsize;
real gamma,dderiv, dderiv1;
// Directional derivative along dir
dderiv = 2.0*dot(nstates,Ygrad,Ydir) + dot(nstates,Bgrad,Bdir,subrot);
// Try to do line minimzations until you succeed or the limit is reached
#define DO_LINMIN_LIMIT 5
int i;
for (i = 0; i < DO_LINMIN_LIMIT; i++){
// Shift Y by stepsize: Y[q] += stepsize*Ydir[q]
scale_accumulate(nstates,stepsize,Ydir,Y);
scale_accumulate(nstates,stepsize,Bdir,B,subrot);
// Get the new directional derivative
// Calculate gradient (and precond gradient if needed)
// NOTE: We need to recalculate U,V,C,n,and d because we've shifted Y.
// Not done in calc_grad_pgrad! (and shouldn't be, since it's redundant
// for minin_loop) -- mhe 3/8/02
calc_UVCn_d(e);
calc_grad_pgrad(Ygrad,NULL,Bgrad,NULL,0,e);
dderiv1 = 2.0*dot(nstates,Ygrad,Ydir) + dot(nstates,Bgrad,Bdir,subrot);
// Shift back
scale_accumulate(nstates,-stepsize,Ydir,Y);
scale_accumulate(nstates,-stepsize,Bdir,B,subrot);
// If dderiv is of a different sign then dderiv, we've gone too far:
// "up the other side of the parabola".
// if((dderiv < 0.) && (dderiv1 > 0.))
// {
// stepsize *= 0.5;
// continue;
// }
// // If dderiv1 is larger than dderiv, or to close, take a bigger step
// // !!!!!!!!!!!!!!!!!!!!! HACK WARNING !!!!!!!!!!!!!!!!!!!!!!
// // we should replace the == with a more inteligent check for real numbers
// // !!!!!!!!!!!!!!!!!!!!! HACK WARNING !!!!!!!!!!!!!!!!!!!!!!
// if(((stepsize < 0.) && (dderiv1 > dderiv) ||
// ((stepsize > 0.) && (dderiv1 < dderiv))) || dderiv1 == dderiv)
// {
// stepsize *= 4.0;
// continue;
// }
// Do a linear fit to dderiv and dderiv1
// to get proposed minimum (gamma)
gamma = (dderiv*stepsize)/(dderiv-dderiv1);
dft_log("dderiv = %25.15le dderiv1 = %25.15le\n",dderiv,dderiv1);
dft_log("stepsize = %14.7le gamma = %14.7le\n",stepsize,gamma);
dft_log_flush();
if(gamma > 0.)
{
gamma = stepsize;
break;
}
// dft_log("dderiv = %25.15le dderiv1 = %25.15le\n",dderiv,dderiv1);
// If the proposed minimum is much larger than the stepsize,
// increase stepsize
if (fabs(gamma/stepsize) > 10.0){
stepsize *= 10.0;
continue;
}
// If much smaller, decrease
// Tairan (10/15/1998): but the stepsize should not be too small.
// otherwise, numerical error will overwhelm the calculation
// Y += stepsize * dir
else if ( fabs(gamma/stepsize)<0.1 &&
fabs(stepsize)>e.cntrl.stepsize_min ){
stepsize *= 0.1;
continue;
}
// Otherwise, it was a good linmin so stop!
else
break;
}
#ifdef DFT_PROFILING
timerOff(9); // Turn off do_linmin timer
#endif // DFT_PROFILING
// see if linmin failed and return 0.0 if it did.
// calling routine should catch this
if ( i == DO_LINMIN_LIMIT )
return 0.0;
// If update_stepsize flag is set, use gamma as the new stepsize
// if it is larger than the stepsize minimum.
// Otherwise, update stepsize only if it has been updated by
// factors of 10 or 0.1
if (e.cntrl.update_stepsize && fabs(gamma)>e.cntrl.stepsize_min){
dft_log("stepsize updated to %e : (min is %e)\n",
gamma, e.cntrl.stepsize_min);
dft_log_flush();
e.cntrl.elec_stepsize = gamma;
}
// Return optimum step found
return gamma;
}
//
// Computes the gradient and possibly the preconditioned gradient
// of the total electronic energy
//
static void
calc_electronic_grad_pgrad(ColumnBundle **Ygrad, ColumnBundle **pYgrad,
Matrix **Bgrad, Matrix **pBgrad,
int precond, Everything &e)
{
if (precond){
dft_log("Computing gradient, preconditioned gradient, and Hsub\n");
dft_log_flush();
calc_elecgrad_pgrad_and_Hsub(Ygrad,pYgrad,Bgrad,pBgrad,e);
}
else {
dft_log("Computing gradient and Hsub\n");
dft_log_flush();
calc_elecgrad_and_Hsub(Ygrad,Bgrad,e);
}
}
//
// Computes the gradient and possibly the preconditioned gradient
// of the band structure energy
//
static void
calc_band_struct_grad_pgrad(ColumnBundle **Ygrad, ColumnBundle **pYgrad,
Matrix **Bgrad, Matrix **pBgrad,
int precond, Everything &e)
{
if (precond)
{
dft_log("Computing gradient and preconditioned gradient\n");
dft_log_flush();
calc_elecgrad_pgrad_band_struct(Ygrad,pYgrad,e);
}
else
{
dft_log("Computing gradient and Hsub\n");
dft_log_flush();
calc_elecgrad_band_struct(Ygrad,e);
}
}
//
// Print gradient length information
//
static void
print_gradient(real abs2g,real abs2Yg,int subrot,real abs2Bg)
{
dft_log("|grad| = %le |Ygrad| = %le",sqrt(abs2g),sqrt(abs2Yg));
if (subrot)
dft_log(" |Bgrad| = %le\n",sqrt(abs2Bg));
dft_log("\n");
dft_log_flush();
}
//
// The main minimization loop: this routine is what does the actual
// work. It is slightly convoluted looking because it must do all
// possibilities of flags: subspace rotations, preconditioning,
// doing line-minimizations, and doing a conjugate gradients algorithm.
//
// It assumes that the space in Ygrad, pYgrad, and Ydir (and their B
// counterparts if subspace rotation is on) is already allocated if
// required by the current set of flags.
//
// It is passed 3 functions: how to calculate the energy which
// is to be minimized, how to print the energy, and how to compute
// its gradient (with possible preconditioning).
//
static void
minim_loop(int niter,
int precond, int subrot,
int do_linmin, int cg_algorithm,
real (*calc_energy)(Everything &e,int inter),
void (*print_energy)(Everything &e,int inter),
void (*calc_grad_pgrad)(ColumnBundle **Ygrad,ColumnBundle **pYgrad,
Matrix **Bgrad,Matrix **pBgrad,
int precond,Everything &e),
ColumnBundle **Ygrad, ColumnBundle **pYgrad, ColumnBundle **Ydir,
Matrix **Bgrad, Matrix **pBgrad, Matrix **Bdir,
Everything &e)
{
// Some local variables, hopefully the names are pretty self
// explanitory. "gamma" is how much we will move along the current
// serach direction. alpha is the amount of the old direction mixed
// in with the current gradient when computing a search direction in
// a conjugate gradient algorithm.
real abs2g, abs2pg, gdotpg, golddotpgold=0.0;
real abs2Yg, abs2Bg;
real gdotdir=0.0, golddotdir=0.0, linminfact=0.0, gamma=0.0;
real alpha;
Control &cntrl = e.cntrl;
// If we are not preconditioning, then point the preconditioned
// gradients to the real gradients (precond by identity). This
// will make the code simpler below.
if (!precond) {
pYgrad = Ygrad;
pBgrad = Bgrad;
}
// If we are not doing a conjugate gradient algorithm, then
// our search direction is along the preconditioned gradient
if (!cg_algorithm) {
Ydir = pYgrad;
Bdir = pBgrad;
}
// Get the current elec variables
int nstates = e.elecinfo.nstates;
ColumnBundle **Y = e.elecvars.Y;
Matrix **B = e.elecvars.B;
// HACK ALERT!!!
// Initial value for global_abs2Yg.
global_abs2Yg = 0.001;
// Here we go!
real Eprev, E = 0.;
real ddotg,magd2,dcosine;
for (int iter=0; iter < niter; iter++) {
// Calculate the current energy
Eprev = E;
E = calc_energy(e,iter);
// Monitor performance of line minimizations (if we are doing them)
if (do_linmin && iter>0)
golddotdir = 2*dot(nstates,Ygrad,Ydir) +
dot(nstates,Bgrad,Bdir,subrot);
// Calculate gradient (and precond gradient if needed)
calc_grad_pgrad(Ygrad,pYgrad,Bgrad,pBgrad,precond,e);
// Square length of gradient, of precond. gradient,
// dot product of gradient precond gradient
abs2Yg = 4.0*abs2(nstates,Ygrad);
abs2Bg = abs2(nstates,Bgrad,subrot);
// HACK ALERT!!!
// Store the magnitude of the gradient in a global variable, so the
// wavelet invL routine can calculate the number of digits to converge
// to.
global_abs2Yg = abs2Yg;
// END OF HACK!!!
abs2g = abs2Yg + abs2Bg;
if (precond) {
abs2pg = 4.0*abs2(nstates,pYgrad) + abs2(nstates,pBgrad,subrot);
gdotpg = 4.0*dot(nstates,Ygrad,pYgrad)+ dot(nstates,Bgrad,pBgrad,subrot);
} else
abs2pg = gdotpg = abs2g;
// DEBUGGING CODE: output the cosine of the angle between Ydir
// and Ygrad. Should be close to zero. -- mhe 02/28/02
if(iter > 0)
{
ddotg = dot(nstates,Ydir,Ygrad);
magd2 = dot(nstates,Ydir,Ydir);
//magg2 = dot(nstates,Ygrad,Ygrad); --> use abs2Yg.
dcosine = ddotg/sqrt(magd2*abs2Yg);
dft_log("\nIteration %d: cos(angle between Ydir and Ygrad) = %1.12lf\n",
iter,dcosine);
}
// END DEBUGGING CODE
// If we are doing line minimizations, then figure
// out how good the previous line min was. We can only
// do this after we've done an iteration.
if (do_linmin && iter>0) {
gdotdir = 2*dot(nstates,Ygrad,Ydir) + dot(nstates,Bgrad,Bdir,subrot);
linminfact = gdotdir / golddotdir;
dft_log("\nlinmin = %le\n\n",linminfact);
dft_log_flush();
}
// Print current status: energy, band structure, and gradient info.
print_energy(e,iter);
print_gradient(abs2g,abs2Yg,subrot,abs2Bg);
// Are we done?
// If this is the last iteration, we can skip all the work below...
if (iter==niter-1) {
// If energy diffs are not converged to set limit, say something
if (!e.cntrl.if_e_converged(e.energies.Etot))
e.cntrl.print(iter);
// Reset convergence tester and get out
e.cntrl.reset_converge();
break;
}
// If energy differences are converged to our tolerance
else if (e.cntrl.if_e_converged(e.energies.Etot)) {
// Reset the convergence tester and get out
e.cntrl.reset_converge();
break;
}
// if we need to, check to see if force differences have converged
else if (e.cntrl.force_diff_check != 0 &&
((iter+1)%e.cntrl.force_diff_check) == 0 &&
e.cntrl.if_force_diff_converged(e)) {
// Reset convergence tester and get out
e.cntrl.reset_converge();
break;
}
// For conjugate gradient methods,
// alpha is how much of old direction gets mixed
// with current gradient to find the new search direction
if (cg_algorithm) {
if (iter>0 && fabs(linminfact) < 0.05)
alpha = gdotpg / golddotpgold;
else
alpha = 0.0;
if(iter>0)
dft_log("alpha = %lg, factor = %lg\n",alpha,
(gdotpg/golddotpgold));
else
dft_log("alpha = %lf\n",alpha);
dft_log_flush();
}
// If we are doing congjugate gradients, then we must
// compute the search direction
// Ydir <- alpha*Ydir + 2.0*[p]Ygrad
if (cg_algorithm) {
scaled_sum(nstates,2.0,pYgrad,alpha,Ydir,Ydir);
scaled_sum(nstates,1.0,pBgrad,alpha,Bdir,Bdir,subrot);
}
// Do linmin if needed to find gamma, which is how much
// along the search direction we are supposed to move
if (do_linmin){
switch(cntrl.linmin_method){
case QUAD:
gamma = linmin(Ydir,Ygrad,subrot,Bdir,Bgrad,calc_energy,e);
break;
case LIN:
gamma = linmin_derivs_only(Ydir,Ygrad,subrot,Bdir,Bgrad,calc_grad_pgrad,e);
break;
default:
die("Invalid linmin_method\n");
}
// check to see if linmin failed
if( gamma == 0.0 ){
dft_log("\n***** WARNING *******\nLINMIN FAILED!!!\n");
dft_log("EXITING MINIMIZATION LOOP\n\n");
dft_log_flush();
break;
}
}
// Without a linmin, then if we did well the previous
// iteration, increase gamma; or decrease if we did badly.
else if (iter>0)
{
if (E < Eprev)
gamma *= 1.05;
else
gamma *= 0.3;
}
// Otherwise, use the stepsize passed in
else
gamma = e.cntrl.elec_stepsize;
// Shift along by gamma
dft_log("Shifting by gamma = %le\n",gamma);
dft_log_flush();
scale_accumulate(nstates,gamma,Ydir,Y);
scale_accumulate(nstates,gamma,Bdir,B,subrot);
// We need these quantities for the next iteration
golddotpgold = gdotpg;
} // iter for loop
}
//
// The routine called from the outside world to do the minimization.
// It figures out the various minimization options that are asked for,
// allocates memory for doing the minimization, and then calls
// the minimization loop function above to do the real work.
//
void
minimize_elec(Everything &e)
{
// What kind of algorithm are we running?
int subrot = FALSE;
int precond = FALSE;
int do_linmin = FALSE;
int cg_algorithm = FALSE;
if (e.elecinfo.subspace_rotation)
subrot = TRUE;
switch(e.cntrl.algorithm){
case DFT_CG:
cg_algorithm = do_linmin = TRUE;
break;
case DFT_PCG:
cg_algorithm = do_linmin = precond = TRUE;
break;
case DFT_SD:
do_linmin = TRUE;
break;
case DFT_PSD:
do_linmin = precond = TRUE;
break;
case DFT_EOM:
break;
case DFT_PEOM:
precond = TRUE;
break;
// We should not get here, but just in case...
default:
die("Unknown minimization algorithm in minimize_elec()!\n");
break;
}
// number of iterations
int niter = e.cntrl.max_elecs_steps;
// Print a banner saying who we are and what we are doing
dft_log("\n----- minimize_elec() -----\n");
time_t timenow = time(0);
dft_log("Called on %s",ctime(&timenow));
dft_log("Starting %d iterations of minimization\n",niter);
if (subrot)
dft_log("* Subspace rotation on\n");
if (precond)
dft_log("* Preconditioning on\n");
if (do_linmin)
dft_log("* Performing line minimizations\n");
if (cg_algorithm)
dft_log("* Using a conjugate gradient algorithm\n");
dft_log("Initial electronic stepsize = %lg\n",e.cntrl.elec_stepsize);
if (e.elecinfo.calc_fillings)
dft_log("Calculating Fermi fillings every %d iterations with kT = %lg\n",
e.elecinfo.niter_recalc_fillings,e.elecinfo.kT);
dft_log("\n");
dft_log_flush();
//
// Allocate memory for what we need to do depending on flags:
//
dft_log("Allocating workspace for minimization\n\n");
dft_log_flush();
// First set everything to NULL
ColumnBundle **Ygrad = NULL;
ColumnBundle **pYgrad = NULL;
ColumnBundle **Ydir = NULL;
Matrix **Bgrad = NULL;
Matrix **pBgrad = NULL;
Matrix **Bdir = NULL;
// Now, get some basic information for the allocation
int nstates = e.elecinfo.nstates;
int nbands = e.elecinfo.nbands;
BlochState *states = e.elecvars.states;
//ipd - no need Basis *basis = e.basis;
// Allocate space for gradients,
Ygrad = alloc_ColumnBundle_array(nstates,states);
if (subrot)
Bgrad = alloc_Matrix_array(nstates,nbands,nbands);
// preconditioned gradients, and
if (precond) {
pYgrad = alloc_ColumnBundle_array(nstates,states);
if (subrot)
pBgrad = alloc_Matrix_array(nstates,nbands,nbands);
}
// search directions (if doing conjugate gradients)
if (cg_algorithm) {
Ydir = alloc_ColumnBundle_array(nstates,states);
if (subrot)
Bdir = alloc_Matrix_array(nstates,nbands,nbands);
}
// Call the minimizer loop to do the hard work
if (e.cntrl.electronic_minimization_flag == TRUE)
minim_loop(niter,precond,subrot,do_linmin,cg_algorithm,
calc_electronic_energy,
print_electronic_energy,
calc_electronic_grad_pgrad,
Ygrad,pYgrad,Ydir,
Bgrad,pBgrad,Bdir,
e);
else if (e.cntrl.band_structure_flag == TRUE)
minim_loop(niter,precond,subrot,do_linmin,cg_algorithm,
calc_band_struct_energy,
print_band_struct_energy,
calc_band_struct_grad_pgrad,
Ygrad,pYgrad,Ydir,
Bgrad,pBgrad,Bdir,
e);
else
die("\nminimize_elec given an unknown minimization method!\n\n");
// Free memory
free_ColumnBundle_array(nstates,Ygrad);
if (subrot)
free_Matrix_array(nstates,Bgrad);
if (precond){
free_ColumnBundle_array(nstates,pYgrad);
if (subrot)
free_Matrix_array(nstates,pBgrad);
}
if (cg_algorithm) {
free_ColumnBundle_array(nstates,Ydir);
if (subrot)
free_Matrix_array(nstates,Bdir);
}
}
syntax highlighted by Code2HTML, v. 0.9.1