ergo
|
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) */