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