ergo
slr.h
Go to the documentation of this file.
1 /* Ergo, version 3.3, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2013 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek.
4  *
5  * This program is free software: you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation, either version 3 of the License, or
8  * (at your option) any later version.
9  *
10  * This program is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  * GNU General Public License for more details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with this program. If not, see <http://www.gnu.org/licenses/>.
17  *
18  * Primary academic reference:
19  * Kohn−Sham Density Functional Theory Electronic Structure Calculations
20  * with Linearly Scaling Computational Time and Memory Usage,
21  * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek,
22  * J. Chem. Theory Comput. 7, 340 (2011),
23  * <http://dx.doi.org/10.1021/ct100611z>
24  *
25  * For further information about Ergo, see <http://www.ergoscf.org>.
26  */
27 
28 #if !defined(SLR_HEADER)
29 #define SLR_HEADER
30 /* Copyright(c) Pawel Salek 2006. */
31 
32 #include <unistd.h>
33 
34 #include "realtype.h"
35 
36 namespace LR {
37 class VarVector;
40 template<bool MultByS,bool SwapXY>
42  public:
43  const VarVector& vec;
45  explicit VarVectorProxyOp(const VarVector& a, ergo_real s=1.0) : vec(a), scalar(s) {}
46 };
47 
48 
52 class VarVector {
54  public:
55  char *fName;
57  int nvar;
58  unsigned onDisk:1;
59  unsigned inMemory:1;
61  VarVector(const VarVector& a) : fName(NULL), nvar(a.nvar),
62  onDisk(0), inMemory(1) {
63  if(nvar) {
64  v = new ergo_real[nvar*2];
65  memcpy(v, a.v, 2*nvar*sizeof(ergo_real));
66  } else v = NULL;
67  }
68 
69  VarVector() : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {}
70 
72  VarVector(int nbast, int nocc, const ergo_real *full_mat)
73  : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {
74  setFromFull(nbast, nocc, full_mat);
75  }
76 
78  if(v) delete v;
79  if(fName) {
80  unlink(fName);
81  delete fName;
82  }
83 }
84  ergo_real* x() const { return v; }
85  ergo_real* y() const { return v + nvar; }
86  void symorth(void);
87  void setSize(int n) {
88  if(nvar != n) {
89  if(v) delete v; v = new ergo_real[2*n];
90  nvar = n;
91  onDisk = 0;
92  }
93  }
94  int size() const { return nvar; }
95  void print(const char *comment = NULL) {
96  if(comment) puts(comment);
97  for(int i=0; i<nvar*2; i++) printf("%15.10f\n", (double)v[i]);
98  }
99  void setFromFull(int nbast, int nocc, const ergo_real *full_mat);
100  void setFull(int nbast, int nocc, ergo_real *full_mat) const;
101  const ergo_real& operator[](int i) const { return v[i]; }
102  ergo_real& operator[](int i) { return v[i]; }
104  for(int i=0; i<2*nvar; i++) v[i] = scalar;
105  onDisk = 0;
106  return *this;
107  };
109  if(this == &b) return *this;
110  if(nvar != b.nvar) {
111  nvar = b.nvar;
112  if(v) delete v;
113  v = new ergo_real[2*nvar];
114  }
115  memcpy(v, b.v, 2*nvar*sizeof(v[0]));
116  onDisk = 0;
117  return *this;
118  }
119 
120  VarVector&
122  if (&proxy.vec == this)
123  throw "VarVector self-assignment not-implemented";
124  if(nvar != proxy.vec.nvar) {
125  if(v) delete v;
126  nvar = proxy.vec.nvar;
127  v = new ergo_real[2*nvar];
128  }
129 
130  for(int i=0; i<2*nvar; i++)
131  v[i] = proxy.scalar*proxy.vec[i];
132  onDisk = 0;
133  return *this;
134  }
135  VarVector&
137  if (&proxy.vec == this)
138  throw "VarVector self-assignment not-implemented";
139  if(nvar != proxy.vec.nvar) {
140  if(v) delete v;
141  nvar = proxy.vec.nvar;
142  v = new ergo_real[2*nvar];
143  }
144  for(int i=0; i<nvar; i++) {
145  v[i] = proxy.scalar*proxy.vec[i+nvar];
146  v[i+nvar] = proxy.scalar*proxy.vec[i];
147  }
148  onDisk = 0;
149  return *this;
150  }
151 
153  void load(const char* tmpdir);
155  void save(const char* tmpdir);
157  void release(const char* tmpdir);
158 
159  friend ergo_real operator*(const VarVector& a, const VarVector& b);
160  friend ergo_real operator*(const VarVector &a,
162  friend ergo_real operator*(const VarVector &a,
164 };
165 
169 class E2Evaluator {
170  public:
171  virtual bool transform(const ergo_real *dmat, ergo_real *fmat) = 0;
172  virtual ~E2Evaluator() {}
173 };
174 
175 /* ################################################################### */
179  unsigned *ages;
180  unsigned currentAge;
181  int nVecs;
183  bool diskMode;
184  public:
185  explicit VarVectorCollection(int nSize=0) : vecs(NULL), ages(NULL), currentAge(0),
186  nVecs(0), nAllocated(0), diskMode(false) {
187  if(nSize) setSize(nSize);
188  }
190  void setSize(int sz);
191  VarVector& operator[](int i);
192  int size() const { return nVecs; }
193  bool getDiskMode() const { return diskMode; }
194  void setDiskMode(bool x) { diskMode = x; }
196  void release();
198  void releaseAll();
199  static const char *tmpdir;
200 };
201 
202 /* ################################################################### */
205  public:
206  virtual void getOper(ergo_real *result) = 0;
207  virtual ~OneElOperator() {}
208 };
209 
211 class SmallMatrix {
213  int nsize;
214  protected:
215  struct RowProxy {
217  explicit RowProxy(ergo_real *r) : toprow(r) {}
218  ergo_real& operator[](int j) const {
219  //printf(" returning row %i -> %p\n", j, toprow + j);
220  return toprow[j]; }
221  };
222  public:
223  explicit SmallMatrix(int sz) : mat(new ergo_real[sz*sz]), nsize(sz) {}
224  ~SmallMatrix() { delete mat; }
225  const RowProxy operator[](int i) {
226  //printf("Returning column %i -> %p\n", i, mat + i*nsize);
227  return RowProxy(mat + i*nsize); }
228  void expand(int newSize);
229 };
230 
231 
232 /* ################################################################### */
235 class LRSolver {
236  public:
237 
238  LRSolver(int nbast, int nocc,
239  const ergo_real *fock_matrix,
240  const ergo_real *s);
241  virtual ~LRSolver() {/* FIXME: delete esub etc */
242  if(cmo) delete cmo;
243  if(fdiag) delete fdiag;
244  delete xSub;
245  }
246 
253  virtual bool getResidual(VarVectorCollection& residualv) = 0;
254 
258  virtual int getInitialGuess(VarVectorCollection& vecs) = 0;
259 
262  virtual ergo_real getPreconditionerShift(int i) const = 0;
263 
265  virtual void increaseSubspaceLimit(int newSize);
266 
268  bool solve(E2Evaluator& e, bool diskMode = false);
272  protected:
273  static const int MVEC = 200;
282  void getAvMinusFreqSv(ergo_real f, ergo_real *weights,
283  VarVector& r);
284 
289  void projectOnSubspace(const VarVector& full, ergo_real *w)/* const*/;
291  void buildVector(const ergo_real *w, VarVector& full) /* const */;
292 
294  void operToVec(OneElOperator& oper, VarVector& res) const;
295 
299  ergo_real setE2diag(int nbast, int nocc,
300  const ergo_real *fock_matrix,
301  const ergo_real *s);
302 
303  int nbast;
304  int nocc;
306  virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
307  void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const;
308  void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const;
309  private:
318  void load_F_MO(ergo_real *fmat) const;
319  bool lintrans(E2Evaluator& e2, const VarVector& v, VarVector& Av) const;
320 };
321 
322 /* ################################################################### */
325 class SetOfEqSolver : public LRSolver {
328  public:
332  const ergo_real *fock_matrix,
333  const ergo_real *s,
334  ergo_real freq)
335  : LRSolver(nbast, nocc, fock_matrix, s), frequency(freq),
336  rhsSub(new ergo_real[MVEC]) {};
337  void setRHS(OneElOperator& op);
338  virtual ~SetOfEqSolver() { delete rhsSub; }
339  virtual ergo_real getPreconditionerShift(int) const { return frequency; }
340  virtual int getInitialGuess(VarVectorCollection& vecs);
341  virtual bool getResidual(VarVectorCollection& residualv);
343  virtual void increaseSubspaceLimit(int newSize);
344  ergo_real getPolarisability(OneElOperator& oper) /* const */;
345  protected:
347  virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
348  ergo_real multiplyXtimesVec(const VarVector& rhs) /* const */;
350 };
351 
352 
353 /* ################################################################### */
355 class EigenSolver : public LRSolver {
358  int nStates;
361  public:
363  const ergo_real *fock_matrix,
364  const ergo_real *overlap, int n)
365  : LRSolver(nbast, nocc, NULL, NULL),
367  nStates(n), nConverged(0), last_ev(NULL) {
368  ritzVals[0] = setE2diag(nbast, nocc, fock_matrix, overlap);
369  for(int i=1; i<n; i++) ritzVals[i] = ritzVals[0];
370  }
371  virtual ~EigenSolver() {
372  if(last_ev) delete last_ev;
373  delete ritzVals;
374  delete transMoms2;
375  }
376  virtual ergo_real getPreconditionerShift(int i) const {
377  return ritzVals[nConverged+i];
378  }
379  virtual int getInitialGuess(VarVectorCollection& vecs);
380  virtual bool getResidual(VarVectorCollection& residualv);
382  virtual void increaseSubspaceLimit(int newSize);
383 
384  ergo_real getFreq(int i) const { return ritzVals[i]; }
385  void computeMoments(OneElOperator& dipx,
386  OneElOperator& dipy,
387  OneElOperator& dipz);
388  ergo_real getTransitionMoment2(int i) const { return transMoms2[i]; }
389 };
390 
391 } /* End of namespace LR */
392 #endif /* !defined(SLR_HEADER) */
void setFromFull(int nbast, int nocc, const ergo_real *full_mat)
Definition: slr.cc:400
VarVector(const VarVector &a)
Definition: slr.h:61
virtual void increaseSubspaceLimit(int newSize)
expands above the default limit
Definition: slr.cc:731
virtual bool transform(const ergo_real *dmat, ergo_real *fmat)=0
void getAvMinusFreqSv(ergo_real f, ergo_real *weights, VarVector &r)
Computes a vector built of base vectors with specified vectors.
Definition: slr.cc:848
LRSolver(int nbast, int nocc, const ergo_real *fock_matrix, const ergo_real *s)
Initialize the solver by computing the diagonal of the E2 operator as needed for preconditioning.
Definition: slr.cc:385
ergo_real * last_ev
most recent eigenvectors in the reduced space
Definition: slr.h:360
virtual void increaseSubspaceLimit(int newSize)
expands above the default limit
Definition: slr.cc:1169
EigenSolver(int nbast, int nocc, const ergo_real *fock_matrix, const ergo_real *overlap, int n)
Definition: slr.h:362
void setRHS(OneElOperator &op)
initializes the rhs field
Definition: slr.cc:1004
double ergo_real
Definition: realtype.h:53
ergo_real setE2diag(int nbast, int nocc, const ergo_real *fock_matrix, const ergo_real *s)
setE2diag is called by the constructor to fill in the approximation of the E[2] operator diagonal...
Definition: slr.cc:800
const ergo_real & operator[](int i) const
Definition: slr.h:101
int nAllocated
Definition: slr.h:182
virtual void addToSpace(VarVectorCollection &vecs, E2Evaluator &e2)
extends the subspace with v and its transformed vector Av.
Definition: slr.cc:1018
Iterative Set Of Linear Equations solver, extending the generic LRSolver.
Definition: slr.h:325
ergo_real * y() const
returns the Y part.
Definition: slr.h:85
RowProxy(ergo_real *r)
Definition: slr.h:217
ergo_real * mat
Definition: slr.h:212
Linear Response iterative solver using a variant of the Davidson method.
Definition: slr.h:235
void setDiskMode(bool x)
Definition: slr.h:194
ergo_real getFreq(int i) const
Definition: slr.h:384
E2Evaluator interface provides a way to perform a linear transformation of supplied transition densit...
Definition: slr.h:169
virtual int getInitialGuess(VarVectorCollection &vecs)
returns the initial guess for the linear set of equations.
Definition: slr.cc:922
const VarVector & vec
Definition: slr.h:43
virtual ~OneElOperator()
Definition: slr.h:207
unsigned inMemory
valid representation in memory
Definition: slr.h:59
virtual int getInitialGuess(VarVectorCollection &vecs)
generate the starting guess for the HOMO-LUMO excitation by placing one in th the right position...
Definition: slr.cc:1188
VarVector & operator=(const VarVectorProxyOp< false, false > &proxy)
Definition: slr.h:121
ergo_real * rhsSub
RHS vector projected onto subspace.
Definition: slr.h:346
VarVectorProxyOp(const VarVector &a, ergo_real s=1.0)
Definition: slr.h:45
void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const
computes mo := cmo'*ao*cmo
Definition: slr.cc:611
VarVectorCollection Avects
vects and Avects members store the trial vectors and their transformed versions.
Definition: slr.h:313
unsigned currentAge
Definition: slr.h:180
void symorth(void)
Uses symmetric orthogonalization to orthonormalize itself (x y) with the swapped vector (y x)...
Definition: slr.cc:462
ergo_real * ritzVals
recent ritz values in the subspace.
Definition: slr.h:356
template based proxy object that uses bool-valued policies to perform the assignments.
Definition: slr.h:41
VarVectorCollection(int nSize=0)
Definition: slr.h:185
ergo_real & operator[](int i)
Definition: slr.h:102
void load(const char *tmpdir)
Load the object to memory.
Definition: slr.cc:215
Vector of variables parametrising the solution to the linear response equations.
Definition: slr.h:52
void computeMoments(OneElOperator &dipx, OneElOperator &dipy, OneElOperator &dipz)
Definition: slr.cc:1222
VarVector & operator=(ergo_real scalar)
Definition: slr.h:103
int nVecs
Definition: slr.h:181
ergo_real scalar
Definition: slr.h:44
SetOfEqSolver(int nbast, int nocc, const ergo_real *fock_matrix, const ergo_real *s, ergo_real freq)
Creates the set-of-equations solver.
Definition: slr.h:331
virtual ~SetOfEqSolver()
Definition: slr.h:338
void load_F_MO(ergo_real *fmat) const
Definition: slr.cc:621
void expand(int newSize)
increase the dimension of the matrix without losing the data.
Definition: slr.cc:370
virtual ergo_real getPreconditionerShift(int) const
returns the preconditioning shift.
Definition: slr.h:339
Iterative Eigenvalue solver, extending the generic LRSolver.
Definition: slr.h:355
virtual int getInitialGuess(VarVectorCollection &vecs)=0
Computes the initial vector the subspace is to be seeded with.
virtual void addToSpace(VarVectorCollection &vecs, E2Evaluator &e2)
extends the subspace with v and its transformed vector Av.
Definition: slr.cc:637
~SmallMatrix()
Definition: slr.h:224
a collection of vectors, usually handled at once.
Definition: slr.h:177
~VarVectorCollection()
Definition: slr.cc:298
ergo_real * cmo
the MO coefficients.
Definition: slr.h:316
VarVector & operator=(const VarVectorProxyOp< false, true > &proxy)
Definition: slr.h:136
void release()
Make sure there is space for one vector.
Definition: slr.cc:326
void computeExactE2Diag(E2Evaluator &e2)
Definition: slr.cc:784
unsigned onDisk
valid representation on disk
Definition: slr.h:58
ergo_real xTimesRHS
Definition: slr.h:349
const RowProxy operator[](int i)
Definition: slr.h:225
VarVector(int nbast, int nocc, const ergo_real *full_mat)
Creates a vector from a full matrix.
Definition: slr.h:72
ergo_real * toprow
Definition: slr.h:216
void print(const char *comment=NULL)
Definition: slr.h:95
virtual void getOper(ergo_real *result)=0
ergo_real getPolarisability(OneElOperator &oper)
computes polarizability by contracting the response vector with specified operator ...
Definition: slr.cc:1035
ergo_real frequency
frequency for which the SOE is to be solved.
Definition: slr.h:326
~VarVector()
Definition: slr.h:77
SmallMatrix(int sz)
Definition: slr.h:223
void releaseAll()
Release all vectors from the memory, saving if necessary.
Definition: slr.cc:352
ergo_real * x() const
return the X part of the vector.
Definition: slr.h:84
a class implementing dynamic resized two dimensional arrays.
Definition: slr.h:211
char * fName
Present in secondary storage (i.e.
Definition: slr.h:55
ergo_real * transMoms2
most recent SQUARED transition moments.
Definition: slr.h:357
VarVector rhs
RHS of the SOE.
Definition: slr.h:327
int maxSubspaceSize
current subspace size limit.
Definition: slr.h:271
static const int MVEC
default limit for subspace size
Definition: slr.h:273
Abstract interface to a one electron operator.
Definition: slr.h:204
int nsize
Definition: slr.h:213
void setSize(int n)
Definition: slr.h:87
int size() const
Definition: slr.h:192
static const char * tmpdir
Definition: slr.h:199
void projectOnSubspace(const VarVector &full, ergo_real *w)
Projects vector.
Definition: slr.cc:872
VarVector & operator[](int i)
Definition: slr.cc:315
int nStates
number of excited states to compute
Definition: slr.h:358
ergo_real convThreshold
iterative method convergence threshold
Definition: slr.h:270
void setFull(int nbast, int nocc, ergo_real *full_mat) const
Definition: slr.cc:429
bool getDiskMode() const
Definition: slr.h:193
void save(const char *tmpdir)
Save the object.
Definition: slr.cc:246
int subspaceSize
current subspace size
Definition: slr.h:275
void setSize(int sz)
Definition: slr.cc:305
void buildVector(const ergo_real *w, VarVector &full)
Build full fector from the reduced form.
Definition: slr.cc:885
ergo_real * fdiag
the eigenvalues of the Fock matrix.
Definition: slr.h:314
SmallMatrix sSub
S[2] matrix projected onto subspace.
Definition: slr.h:277
void operToVec(OneElOperator &oper, VarVector &res) const
Transform square operator to the vector form.
Definition: slr.cc:898
int size() const
Definition: slr.h:94
int nConverged
number of already converged eigenvalues
Definition: slr.h:359
virtual ergo_real getPreconditionerShift(int i) const
returns the preconditioning shift.
Definition: slr.h:376
virtual void increaseSubspaceLimit(int newSize)
expands above the default limit
Definition: slr.cc:992
VarVector e2diag
approximation to the diagonal of E2 operator
Definition: slr.h:274
Definition: slr.h:215
bool diskMode
Definition: slr.h:183
VarVector & operator=(const VarVector &b)
Definition: slr.h:108
virtual ergo_real getPreconditionerShift(int i) const =0
returns the preconditioning shift.
SmallMatrix eSub
E[2] matrix projected onto subspace.
Definition: slr.h:276
bool lintrans(E2Evaluator &e2, const VarVector &v, VarVector &Av) const
performs the linear transformation of the vector with E[2] operator.
Definition: slr.cc:703
ergo_real & operator[](int j) const
Definition: slr.h:218
VarVector()
Definition: slr.h:69
int nvar
nvar := nocc*nvirt.
Definition: slr.h:57
void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const
Definition: slr.cc:600
ergo_real * v
vector coefficients
Definition: slr.h:53
ergo_real multiplyXtimesVec(const VarVector &rhs)
multiplies current solution by some vector.
Definition: slr.cc:946
VarVectorCollection vects
base vectors
Definition: slr.h:305
friend ergo_real operator*(const VarVector &a, const VarVector &b)
Definition: slr.cc:97
ergo_real getTransitionMoment2(int i) const
Definition: slr.h:388
int nbast
number of basis functions
Definition: slr.h:303
int nocc
number of occupied orbitals
Definition: slr.h:304
virtual bool getResidual(VarVectorCollection &residualv)
get residual of the eigenvalue problem.
Definition: slr.cc:1054
unsigned * ages
Definition: slr.h:179
virtual ~E2Evaluator()
Definition: slr.h:172
ergo_real * xSub
solution vector projected onto subspace
Definition: slr.h:278
void release(const char *tmpdir)
Releases the memory, saving if necessary.
Definition: slr.cc:282
virtual ~EigenSolver()
Definition: slr.h:371
virtual bool getResidual(VarVectorCollection &residualv)=0
Computes the residual vector.
bool solve(E2Evaluator &e, bool diskMode=false)
Solves the problem defined by the subclass.
Definition: slr.cc:749
virtual ~LRSolver()
Definition: slr.h:241
VarVector * vecs
Definition: slr.h:178
virtual bool getResidual(VarVectorCollection &residualv)
get the residual of the set of linear equations.
Definition: slr.cc:963