Pulay.h

Go to the documentation of this file.

/*-------------------------------------------------------------------
Copyright 2015 Ravishankar Sundararaman

This file is part of JDFTx.

JDFTx 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 3 of the License, or
(at your option) any later version.

JDFTx 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 JDFTx.  If not, see <http://www.gnu.org/licenses/>.
-------------------------------------------------------------------*/

#ifndef JDFTX_CORE_PULAY_H
#define JDFTX_CORE_PULAY_H

#include <core/PulayParams.h>
#include <core/Util.h>
#include <electronic/matrix.h>
#include <cfloat>

template<typename Variable> class Pulay
{
public:
        Pulay(const PulayParams& pp);
        
        double minimize(double Eprev=+DBL_MAX, std::vector<string> extraNames=std::vector<string>(), std::vector<double> extraThresh=std::vector<double>());
        
        void loadState(const char* filename); 
        void saveState(const char* filename) const; 

        virtual double sync(double x) const { return x; }
        
protected:
        //----- Interface specification -----

        virtual double cycle(double dEprev, std::vector<double>& extraValues)=0;
        
        virtual void report(int iter) {} 
        virtual void axpy(double alpha, const Variable& X, Variable& Y) const=0; 
        virtual double dot(const Variable& X, const Variable& Y) const=0; 
        virtual size_t variableSize() const=0; 
        virtual void readVariable(Variable&, FILE*) const=0; 
        virtual void writeVariable(const Variable&, FILE*) const=0; 
        virtual Variable getVariable() const=0; 
        virtual void setVariable(const Variable&)=0; 
        virtual Variable precondition(const Variable&) const=0; 
        virtual Variable applyMetric(const Variable&) const=0; 

private:
        const PulayParams& pp; 
        std::vector<Variable> pastVariables; 
        std::vector<Variable> pastResiduals; 
        matrix overlap; 
};


//---------------------- Implementation ----------------------------

#include <core/Minimize_linmin.h>
#include <memory>

//Norm convergence check (eigenvalue-difference or residual)
//Make sure value is within tolerance for nCheck consecutive cycles
class NormCheck
{       unsigned nCheck; double threshold;
        std::deque<bool> history;
public:
        NormCheck(unsigned nCheck, double threshold) : nCheck(nCheck), threshold(fabs(threshold)) {}
        bool checkConvergence(double norm)
        {       history.push_back(fabs(norm)<threshold);
                if(history.size()==nCheck+1) history.pop_front(); //discard old unneeded elements 
                if(history.size()==nCheck)
                {       for(bool converged: history)
                                if(!converged)
                                        return false;
                        return true;
                }
                else return false;
        }
};

template<typename Variable> Pulay<Variable>::Pulay(const PulayParams& pp)
: pp(pp), overlap(pp.history, pp.history)
{
}

template<typename Variable> double Pulay<Variable>::minimize(double Eprev, std::vector<string> extraNames, std::vector<double> extraThresh)
{
        double E = sync(Eprev); Eprev = 0.;
        double dE = E-Eprev;
        assert(extraNames.size()==extraThresh.size());
        
        //Initialize convergence checkers:
        EdiffCheck ediffCheck(2, pp.energyDiffThreshold); ediffCheck.checkConvergence(E); //store the initial energy in the check's history
        NormCheck resCheck(2, pp.residualThreshold);
        std::vector<std::shared_ptr<NormCheck> > extraCheck(extraNames.size());
        for(size_t iExtra=0; iExtra<extraNames.size(); iExtra++)
                extraCheck[iExtra] = std::make_shared<NormCheck>(2, extraThresh[iExtra]);

        for(int iter=0; iter<pp.nIterations; iter++)
        {
                //If history is full, remove oldest member
                assert(pastResiduals.size() == pastVariables.size());
                if((int)pastResiduals.size() >= pp.history)
                {       size_t ndim = pastResiduals.size();
                        if(ndim>1) overlap.set(0,ndim-1, 0,ndim-1, overlap(1,ndim, 1,ndim));
                        pastVariables.erase(pastVariables.begin());
                        pastResiduals.erase(pastResiduals.begin());
                }
                
                //Cache the old energy and variables
                Eprev = E;
                pastVariables.push_back(getVariable());

                //Perform cycle:
                std::vector<double> extraValues(extraThresh.size());
                E = sync(cycle(dE, extraValues));
                dE = E - Eprev;
                for(auto& v: extraValues) v = sync(v);
                        
                //Calculate and cache residual:
                double residualNorm = 0.;
                {       Variable residual = getVariable(); axpy(-1., pastVariables.back(), residual);
                        pastResiduals.push_back(residual);
                        residualNorm = sync(sqrt(dot(residual,residual)));
                }
                
                //Print energy and convergence parameters:
                fprintf(pp.fpLog, "%sCycle: %2i   %s: ", pp.linePrefix, iter, pp.energyLabel);
                fprintf(pp.fpLog, pp.energyFormat, E);
                fprintf(pp.fpLog, "   d%s: %+.3e", pp.energyLabel, dE);
                fprintf(pp.fpLog, "   |Residual|: %.3e", residualNorm);
                for(size_t iExtra=0; iExtra<extraNames.size(); iExtra++)
                        fprintf(pp.fpLog, "   |%s|: %.3e", extraNames[iExtra].c_str(), extraValues[iExtra]);
                fprintf(pp.fpLog, "\n"); fflush(pp.fpLog);
                
                //Optional reporting:
                report(iter);
                
                //Check for convergence and update variable:
                bool converged = false;
                if(!converged && ediffCheck.checkConvergence(E))
                {       fprintf(pp.fpLog, "%sConverged (|Delta E|<%le for 2 iters).\n\n", pp.linePrefix, pp.energyDiffThreshold);
                        converged = true;
                }
                if(!converged && resCheck.checkConvergence(residualNorm))
                {       fprintf(pp.fpLog, "%sConverged (|Residual|<%le for 2 iters).\n\n", pp.linePrefix, pp.residualThreshold);
                        converged = true;
                }
                if(!converged && extraNames.size())
                        for(size_t iExtra=0; iExtra<extraNames.size(); iExtra++)
                                if(extraCheck[iExtra]->checkConvergence(extraValues[iExtra]))
                                {       fprintf(pp.fpLog, "%sConverged (|%s|<%le for 2 iters).\n\n", pp.linePrefix, extraNames[iExtra].c_str(), extraThresh[iExtra]);
                                        converged = true;
                                        break;
                                }
                fflush(pp.fpLog);
                if(converged || killFlag) break; //converged or manually interrupted
                
                //---- DIIS/Pulay mixing -----
                        
                //Update the overlap matrix
                size_t ndim = pastResiduals.size();
                Variable MlastResidual = applyMetric(pastResiduals.back());
                for(size_t j=0; j<ndim; j++)
                {       double thisOverlap = dot(pastResiduals[j], MlastResidual);
                        overlap.set(j, ndim-1, thisOverlap);
                        overlap.set(ndim-1, j, thisOverlap);
                }
                
                //Invert the residual overlap matrix to get the minimum of residual
                matrix cOverlap(ndim+1, ndim+1); //Add row and column to enforce normalization constraint
                cOverlap.set(0, ndim, 0, ndim, overlap(0, ndim, 0, ndim));
                for(size_t j=0; j<ndim; j++)
                {       cOverlap.set(j, ndim, 1);
                        cOverlap.set(ndim, j, 1);
                }
                cOverlap.set(ndim, ndim, 0);
                matrix cOverlap_inv = inv(cOverlap);
                
                //Update variable:
                complex* coefs = cOverlap_inv.data();
                Variable v;
                for(size_t j=0; j<ndim; j++)
                {       double alpha = coefs[cOverlap_inv.index(j, ndim)].real();
                        axpy(alpha, pastVariables[j], v);
                        axpy(alpha, precondition(pastResiduals[j]), v);
                }
                setVariable(v);
        }
        return E;
}

template<typename Variable> void Pulay<Variable>::loadState(const char* filename)
{
        size_t nBytesCycle = 2 * variableSize(); //number of bytes per history entry
        size_t nBytesFile = fileSize(filename);
        size_t ndim = nBytesFile / nBytesCycle;
        size_t dimOffset = 0;
        if(int(ndim) > pp.history)
        {       dimOffset = ndim - pp.history;
                ndim = pp.history;
        }
        if(nBytesFile % nBytesCycle != 0)
                die("Pulay history file '%s' does not contain an integral multiple of the mixed variables and residuals.\n", filename);
        fprintf(pp.fpLog, "%sReading %lu past variables and residuals from '%s' ... ", pp.linePrefix, ndim, filename); logFlush();
        pastVariables.resize(ndim);
        pastResiduals.resize(ndim);
        FILE* fp = fopen(filename, "r");
        if(dimOffset) fseek(fp, dimOffset*nBytesCycle, SEEK_SET);
        for(size_t idim=0; idim<ndim; idim++)
        {       readVariable(pastVariables[idim], fp);
                readVariable(pastResiduals[idim], fp);
        }
        fclose(fp);
        fprintf(pp.fpLog, "done.\n"); fflush(pp.fpLog);
        //Compute overlaps of loaded history:
        for(size_t i=0; i<ndim; i++)
        {       Variable Mresidual_i = applyMetric(pastResiduals[i]);
                for(size_t j=0; j<=i; j++)
                {       double thisOverlap = dot(pastResiduals[j], Mresidual_i);
                        overlap.set(i,j, thisOverlap);
                        overlap.set(j,i, thisOverlap);
                }
        }
}

template<typename Variable> void Pulay<Variable>::saveState(const char* filename) const
{
        if(mpiUtil->isHead())
        {       FILE* fp = fopen(filename, "w");
                for(size_t idim=0; idim<pastVariables.size(); idim++)
                {       writeVariable(pastVariables[idim], fp);
                        writeVariable(pastResiduals[idim], fp);
                }
                fclose(fp);
        }
}


#endif //JDFTX_CORE_PULAY_H