slr.h

Go to the documentation of this file.
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure
00002  * calculations.
00003  * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
00004  * 
00005  * This program is free software: you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation, either version 3 of the License, or
00008  * (at your option) any later version.
00009  * 
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  * 
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
00017  * 
00018  * Primary academic reference:
00019  * Kohn−Sham Density Functional Theory Electronic Structure Calculations 
00020  * with Linearly Scaling Computational Time and Memory Usage,
00021  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
00022  * J. Chem. Theory Comput. 7, 340 (2011),
00023  * <http://dx.doi.org/10.1021/ct100611z>
00024  * 
00025  * For further information about Ergo, see <http://www.ergoscf.org>.
00026  */
00027 
00028 #if !defined(SLR_HEADER)
00029 #define SLR_HEADER
00030 /* Copyright(c) Pawel Salek 2006. */
00031 
00032 #include <unistd.h>
00033 
00034 #include "realtype.h"
00035 
00036 namespace LR {
00037 class VarVector;
00040 template<bool MultByS,bool SwapXY>
00041   class VarVectorProxyOp {
00042  public:
00043   const VarVector& vec;
00044   ergo_real scalar;
00045   explicit VarVectorProxyOp(const VarVector& a, ergo_real s=1.0) : vec(a), scalar(s) {}
00046 };
00047 
00048 
00052 class VarVector {
00053   ergo_real *v;  
00054  public:
00055   char *fName;    
00057   int  nvar;      
00058   unsigned onDisk:1;   
00059   unsigned inMemory:1; 
00061   VarVector(const VarVector& a) : fName(NULL), nvar(a.nvar),
00062     onDisk(0), inMemory(1) {
00063     if(nvar) {
00064       v = new ergo_real[nvar*2];
00065       memcpy(v, a.v, 2*nvar*sizeof(ergo_real));
00066     } else v = NULL;
00067   }
00068 
00069   VarVector() : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {}
00070 
00072   VarVector(int nbast, int nocc, const ergo_real *full_mat)
00073     : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1)  {
00074     setFromFull(nbast, nocc, full_mat);
00075   }
00076 
00077   ~VarVector() {
00078     if(v) delete v; 
00079     if(fName) {
00080       unlink(fName);
00081       delete fName;
00082     }
00083 }
00084   ergo_real* x() const { return v; } 
00085   ergo_real* y() const { return v + nvar; } 
00086   void symorth(void);
00087   void setSize(int n) { 
00088     if(nvar != n) {
00089       if(v) delete v; v = new ergo_real[2*n];
00090       nvar = n;
00091       onDisk = 0;
00092     }
00093   }
00094   int size() const { return nvar; }
00095   void print(const char *comment = NULL) {
00096     if(comment) puts(comment);
00097     for(int i=0; i<nvar*2; i++) printf("%15.10f\n", (double)v[i]);
00098   }
00099   void setFromFull(int nbast, int nocc, const ergo_real *full_mat);
00100   void setFull(int nbast, int nocc, ergo_real *full_mat) const;
00101   const ergo_real& operator[](int i) const { return v[i]; }
00102   ergo_real& operator[](int i)             { return v[i]; }
00103   VarVector& operator=(ergo_real scalar) {
00104     for(int i=0; i<2*nvar; i++) v[i] = scalar;
00105     onDisk = 0;
00106     return *this;
00107   };
00108   VarVector& operator=(const VarVector& b) {
00109     if(this == &b) return *this;
00110     if(nvar != b.nvar) { 
00111       nvar = b.nvar; 
00112       if(v) delete v;
00113       v = new ergo_real[2*nvar];
00114     }
00115     memcpy(v, b.v, 2*nvar*sizeof(v[0]));
00116     onDisk = 0;
00117     return *this;
00118   }
00119 
00120     VarVector&
00121     operator=(const VarVectorProxyOp<false, false >& proxy) {
00122       if (&proxy.vec == this)
00123         throw "VarVector self-assignment not-implemented";
00124       if(nvar != proxy.vec.nvar) {
00125         if(v) delete v; 
00126         nvar = proxy.vec.nvar;
00127         v = new ergo_real[2*nvar];
00128       }
00129 
00130       for(int i=0; i<2*nvar; i++)
00131         v[i]      = proxy.scalar*proxy.vec[i];
00132       onDisk = 0;
00133       return *this;
00134     }
00135     VarVector&
00136     operator=(const VarVectorProxyOp<false, true >& proxy) {
00137       if (&proxy.vec == this)
00138         throw "VarVector self-assignment not-implemented";
00139       if(nvar != proxy.vec.nvar) {
00140         if(v) delete v; 
00141         nvar = proxy.vec.nvar;
00142         v = new ergo_real[2*nvar];
00143       }
00144       for(int i=0; i<nvar; i++) {
00145         v[i]      = proxy.scalar*proxy.vec[i+nvar];
00146         v[i+nvar] = proxy.scalar*proxy.vec[i];
00147       }
00148       onDisk = 0;
00149       return *this;
00150     }
00151     
00153     void load(const char* tmpdir);
00155     void save(const char* tmpdir);
00157     void release(const char* tmpdir);
00158 
00159     friend ergo_real operator*(const VarVector& a, const VarVector& b);
00160     friend ergo_real operator*(const VarVector &a,
00161                                const VarVectorProxyOp<false,false>& b);
00162     friend ergo_real operator*(const VarVector &a,
00163                                const VarVectorProxyOp<true,false>& b);
00164 };
00165 
00169 class E2Evaluator {
00170  public:
00171   virtual bool transform(const ergo_real *dmat, ergo_real *fmat) = 0;
00172   virtual ~E2Evaluator() {}
00173 };
00174 
00175 /* ################################################################### */
00177 class VarVectorCollection {
00178   VarVector *vecs;
00179   unsigned *ages;
00180   unsigned currentAge;
00181   int nVecs;
00182   int nAllocated;
00183   bool diskMode;
00184  public:
00185   explicit VarVectorCollection(int nSize=0) : vecs(NULL), ages(NULL), currentAge(0),
00186     nVecs(0), nAllocated(0), diskMode(false) {
00187     if(nSize) setSize(nSize);
00188   }
00189   ~VarVectorCollection();
00190   void setSize(int sz);
00191   VarVector& operator[](int i);
00192   int size() const { return nVecs; }
00193   bool getDiskMode() const { return diskMode; }
00194   void setDiskMode(bool x) { diskMode = x; }
00196   void release(); 
00198   void releaseAll();
00199   static const char *tmpdir;
00200 };
00201  
00202 /* ################################################################### */
00204 class OneElOperator {
00205  public:
00206   virtual void getOper(ergo_real *result) = 0;
00207   virtual ~OneElOperator() {}
00208 };
00209 
00211 class SmallMatrix {
00212   ergo_real *mat;
00213   int nsize;
00214  protected:
00215   struct RowProxy {
00216     ergo_real *toprow;
00217     explicit RowProxy(ergo_real *r) : toprow(r) {}
00218     ergo_real& operator[](int j) const {
00219       //printf("  returning row %i -> %p\n", j, toprow + j);
00220       return toprow[j]; }
00221   };
00222  public:
00223   explicit SmallMatrix(int sz) : mat(new ergo_real[sz*sz]), nsize(sz) {}
00224   ~SmallMatrix() { delete mat; }
00225   const RowProxy operator[](int i) {
00226     //printf("Returning column %i -> %p\n", i, mat + i*nsize);
00227     return RowProxy(mat + i*nsize); }
00228   void expand(int newSize);
00229 };
00230 
00231 
00232 /* ################################################################### */
00235 class LRSolver {
00236  public:
00237 
00238   LRSolver(int nbast, int nocc,
00239            const ergo_real *fock_matrix,
00240            const ergo_real *s);
00241   virtual ~LRSolver() {/* FIXME: delete esub etc */
00242     if(cmo)   delete cmo; 
00243     if(fdiag) delete fdiag;
00244     delete xSub;
00245   }
00246 
00253   virtual bool getResidual(VarVectorCollection& residualv) = 0;
00254 
00258   virtual int getInitialGuess(VarVectorCollection& vecs) = 0;
00259 
00262   virtual ergo_real getPreconditionerShift(int i) const = 0;
00263 
00265   virtual void increaseSubspaceLimit(int newSize);
00266 
00268   bool solve(E2Evaluator& e, bool diskMode = false);
00269   void computeExactE2Diag(E2Evaluator& e2);
00270   ergo_real convThreshold;   
00271   int       maxSubspaceSize; 
00272  protected:
00273   static const int MVEC = 200; 
00274   VarVector e2diag; 
00275   int subspaceSize; 
00276   SmallMatrix eSub; 
00277   SmallMatrix sSub; 
00278   ergo_real *xSub;  
00282   void getAvMinusFreqSv(ergo_real f, ergo_real *weights,
00283                         VarVector& r);
00284 
00289   void projectOnSubspace(const VarVector& full, ergo_real *w)/* const*/;
00291   void buildVector(const ergo_real *w, VarVector& full) /* const */;
00292 
00294   void operToVec(OneElOperator& oper, VarVector& res) const;
00295 
00299   ergo_real setE2diag(int nbast, int nocc,
00300                       const ergo_real *fock_matrix,
00301                       const ergo_real *s);
00302   
00303   int nbast; 
00304   int nocc;  
00305   VarVectorCollection vects;  
00306   virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
00307   void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const;
00308   void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const;
00309  private:
00313   VarVectorCollection Avects; 
00314   ergo_real *fdiag;       
00316   ergo_real *cmo;         
00318   void load_F_MO(ergo_real *fmat) const;
00319   bool lintrans(E2Evaluator& e2, const VarVector& v, VarVector& Av) const;
00320 };
00321 
00322 /* ################################################################### */
00325 class SetOfEqSolver : public LRSolver {
00326   ergo_real frequency; 
00327   VarVector rhs; 
00328  public:
00331   SetOfEqSolver(int nbast, int nocc,
00332                 const ergo_real *fock_matrix,
00333                 const ergo_real *s,
00334                 ergo_real freq)
00335     : LRSolver(nbast, nocc, fock_matrix, s), frequency(freq),
00336     rhsSub(new ergo_real[MVEC]) {};
00337   void setRHS(OneElOperator& op);    
00338   virtual ~SetOfEqSolver() { delete rhsSub; }
00339   virtual ergo_real getPreconditionerShift(int) const { return frequency; }
00340   virtual int getInitialGuess(VarVectorCollection& vecs);
00341   virtual bool getResidual(VarVectorCollection& residualv);
00343   virtual void increaseSubspaceLimit(int newSize);
00344   ergo_real getPolarisability(OneElOperator& oper) /* const */;
00345  protected:
00346   ergo_real *rhsSub; 
00347   virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
00348   ergo_real multiplyXtimesVec(const VarVector& rhs) /* const */;
00349   ergo_real xTimesRHS;
00350 };
00351 
00352 
00353 /* ################################################################### */
00355 class EigenSolver : public LRSolver {
00356   ergo_real *ritzVals; 
00357   ergo_real *transMoms2; 
00358   int nStates;    
00359   int nConverged; 
00360   ergo_real *last_ev; 
00361  public:
00362   EigenSolver(int nbast, int nocc,
00363               const ergo_real *fock_matrix,
00364               const ergo_real *overlap, int n)
00365     : LRSolver(nbast, nocc, NULL, NULL), 
00366     ritzVals(new ergo_real[MVEC]), transMoms2(new ergo_real[MVEC]),
00367     nStates(n), nConverged(0), last_ev(NULL) {
00368     ritzVals[0] = setE2diag(nbast, nocc, fock_matrix, overlap);
00369     for(int i=1; i<n; i++) ritzVals[i] = ritzVals[0];
00370   }
00371   virtual ~EigenSolver() {
00372     if(last_ev) delete last_ev;
00373     delete ritzVals;
00374     delete transMoms2;
00375   }
00376   virtual ergo_real getPreconditionerShift(int i) const { 
00377     return ritzVals[nConverged+i];
00378   }
00379   virtual int getInitialGuess(VarVectorCollection& vecs);
00380   virtual bool getResidual(VarVectorCollection& residualv);
00382   virtual void increaseSubspaceLimit(int newSize);
00383 
00384   ergo_real getFreq(int i) const { return ritzVals[i]; }
00385   void computeMoments(OneElOperator& dipx,
00386                       OneElOperator& dipy,
00387                       OneElOperator& dipz);
00388   ergo_real getTransitionMoment2(int i) const { return transMoms2[i]; }
00389 };
00390 
00391 } /* End of namespace LR */
00392 #endif /* !defined(SLR_HEADER) */

Generated on 21 Nov 2012 for ergo by  doxygen 1.6.1