ergo
slr.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
37
38#if !defined(SLR_HEADER)
39#define SLR_HEADER
40/* Copyright(c) Pawel Salek 2006. */
41
42#include <unistd.h>
43
44#include "realtype.h"
45
46namespace LR {
47class VarVector;
50template<bool MultByS,bool SwapXY>
52 public:
53 const VarVector& vec;
55 explicit VarVectorProxyOp(const VarVector& a, ergo_real s=1.0) : vec(a), scalar(s) {}
56};
57
58
62class VarVector {
64 public:
65 char *fName;
67 int nvar;
68 unsigned onDisk:1;
69 unsigned inMemory:1;
70
71 VarVector(const VarVector& a) : fName(NULL), nvar(a.nvar),
72 onDisk(0), inMemory(1) {
73 if(nvar) {
74 v = new ergo_real[nvar*2];
75 memcpy(v, a.v, 2*nvar*sizeof(ergo_real));
76 } else v = NULL;
77 }
78
79 VarVector() : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {}
80
82 VarVector(int nbast, int nocc, const ergo_real *full_mat)
83 : v(NULL), fName(NULL), nvar(0), onDisk(0), inMemory(1) {
84 setFromFull(nbast, nocc, full_mat);
85 }
86
88 if(v) delete v;
89 if(fName) {
90 unlink(fName);
91 delete fName;
92 }
93}
94 ergo_real* x() const { return v; }
95 ergo_real* y() const { return v + nvar; }
96 void symorth(void);
97 void setSize(int n) {
98 if(nvar != n) {
99 if(v)
100 delete v;
101 v = new ergo_real[2*n];
102 nvar = n;
103 onDisk = 0;
104 }
105 }
106 int size() const { return nvar; }
107 void print(const char *comment = NULL) {
108 if(comment) puts(comment);
109 for(int i=0; i<nvar*2; i++) printf("%15.10f\n", (double)v[i]);
110 }
111 void setFromFull(int nbast, int nocc, const ergo_real *full_mat);
112 void setFull(int nbast, int nocc, ergo_real *full_mat) const;
113 const ergo_real& operator[](int i) const { return v[i]; }
114 ergo_real& operator[](int i) { return v[i]; }
116 for(int i=0; i<2*nvar; i++) v[i] = scalar;
117 onDisk = 0;
118 return *this;
119 };
121 if(this == &b) return *this;
122 if(nvar != b.nvar) {
123 nvar = b.nvar;
124 if(v) delete v;
125 v = new ergo_real[2*nvar];
126 }
127 memcpy(v, b.v, 2*nvar*sizeof(v[0]));
128 onDisk = 0;
129 return *this;
130 }
131
132 VarVector&
134 if (&proxy.vec == this)
135 throw "VarVector self-assignment not-implemented";
136 if(nvar != proxy.vec.nvar) {
137 if(v) delete v;
138 nvar = proxy.vec.nvar;
139 v = new ergo_real[2*nvar];
140 }
141
142 for(int i=0; i<2*nvar; i++)
143 v[i] = proxy.scalar*proxy.vec[i];
144 onDisk = 0;
145 return *this;
146 }
147 VarVector&
149 if (&proxy.vec == this)
150 throw "VarVector self-assignment not-implemented";
151 if(nvar != proxy.vec.nvar) {
152 if(v) delete v;
153 nvar = proxy.vec.nvar;
154 v = new ergo_real[2*nvar];
155 }
156 for(int i=0; i<nvar; i++) {
157 v[i] = proxy.scalar*proxy.vec[i+nvar];
158 v[i+nvar] = proxy.scalar*proxy.vec[i];
159 }
160 onDisk = 0;
161 return *this;
162 }
163
165 void load(const char* tmpdir);
167 void save(const char* tmpdir);
169 void release(const char* tmpdir);
170
171 friend ergo_real operator*(const VarVector& a, const VarVector& b);
172 friend ergo_real operator*(const VarVector &a,
174 friend ergo_real operator*(const VarVector &a,
176};
177
182 public:
183 virtual bool transform(const ergo_real *dmat, ergo_real *fmat) = 0;
184 virtual ~E2Evaluator() {}
185};
186
187/* ################################################################### */
191 unsigned *ages;
192 unsigned currentAge;
193 int nVecs;
196 public:
197 explicit VarVectorCollection(int nSize=0) : vecs(NULL), ages(NULL), currentAge(0),
198 nVecs(0), nAllocated(0), diskMode(false) {
199 if(nSize) setSize(nSize);
200 }
202 void setSize(int sz);
203 VarVector& operator[](int i);
204 int size() const { return nVecs; }
205 bool getDiskMode() const { return diskMode; }
206 void setDiskMode(bool x) { diskMode = x; }
208 void release();
210 void releaseAll();
211 static const char *tmpdir;
212};
213
214/* ################################################################### */
217 public:
218 virtual void getOper(ergo_real *result) = 0;
219 virtual ~OneElOperator() {}
220};
221
225 int nsize;
226 protected:
227 struct RowProxy {
229 explicit RowProxy(ergo_real *r) : toprow(r) {}
230 ergo_real& operator[](int j) const {
231 //printf(" returning row %i -> %p\n", j, toprow + j);
232 return toprow[j]; }
233 };
234 public:
235 explicit SmallMatrix(int sz) : mat(new ergo_real[sz*sz]), nsize(sz) {}
236 ~SmallMatrix() { delete [] mat; }
237 const RowProxy operator[](int i) {
238 //printf("Returning column %i -> %p\n", i, mat + i*nsize);
239 return RowProxy(mat + i*nsize); }
240 void expand(int newSize);
241};
242
243
244/* ################################################################### */
247class LRSolver {
248 public:
249
250 LRSolver(int nbast, int nocc,
251 const ergo_real *fock_matrix,
252 const ergo_real *s);
253 virtual ~LRSolver() {/* FIXME: delete esub etc */
254 if(cmo) delete cmo;
255 if(fdiag) delete fdiag;
256 delete [] xSub;
257 }
258
265 virtual bool getResidual(VarVectorCollection& residualv) = 0;
266
270 virtual int getInitialGuess(VarVectorCollection& vecs) = 0;
271
274 virtual ergo_real getPreconditionerShift(int i) const = 0;
275
277 virtual void increaseSubspaceLimit(int newSize);
278
280 bool solve(E2Evaluator& e, bool diskMode = false);
284 protected:
285 static const int MVEC = 200;
291
294 void getAvMinusFreqSv(ergo_real f, ergo_real *weights,
295 VarVector& r);
296
301 void projectOnSubspace(const VarVector& full, ergo_real *w)/* const*/;
303 void buildVector(const ergo_real *w, VarVector& full) /* const */;
304
306 void operToVec(OneElOperator& oper, VarVector& res) const;
307
311 ergo_real setE2diag(int nbast, int nocc,
312 const ergo_real *fock_matrix,
313 const ergo_real *s);
314
315 int nbast;
316 int nocc;
318 virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
319 void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const;
320 void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const;
321 private:
329
330 void load_F_MO(ergo_real *fmat) const;
331 bool lintrans(E2Evaluator& e2, const VarVector& v, VarVector& Av) const;
332};
333
334/* ################################################################### */
337class SetOfEqSolver : public LRSolver {
340 public:
344 const ergo_real *fock_matrix,
345 const ergo_real *s,
346 ergo_real freq)
347 : LRSolver(nbast, nocc, fock_matrix, s), frequency(freq),
348 rhsSub(new ergo_real[MVEC]) {};
349 void setRHS(OneElOperator& op);
350 virtual ~SetOfEqSolver() { delete [] rhsSub; }
351 virtual ergo_real getPreconditionerShift(int) const { return frequency; }
352 virtual int getInitialGuess(VarVectorCollection& vecs);
353 virtual bool getResidual(VarVectorCollection& residualv);
355 virtual void increaseSubspaceLimit(int newSize);
356 ergo_real getPolarisability(OneElOperator& oper) /* const */;
357 protected:
359 virtual void addToSpace(VarVectorCollection& vecs, E2Evaluator& e2);
360 ergo_real multiplyXtimesVec(const VarVector& rhs) /* const */;
362};
363
364
365/* ################################################################### */
367class EigenSolver : public LRSolver {
373 public:
375 const ergo_real *fock_matrix,
376 const ergo_real *overlap, int n)
377 : LRSolver(nbast, nocc, NULL, NULL),
379 nStates(n), nConverged(0), last_ev(NULL) {
380 ritzVals[0] = setE2diag(nbast, nocc, fock_matrix, overlap);
381 for(int i=1; i<n; i++) ritzVals[i] = ritzVals[0];
382 }
383 virtual ~EigenSolver() {
384 if(last_ev) delete last_ev;
385 delete [] ritzVals;
386 delete [] transMoms2;
387 }
388 virtual ergo_real getPreconditionerShift(int i) const {
389 return ritzVals[nConverged+i];
390 }
391 virtual int getInitialGuess(VarVectorCollection& vecs);
392 virtual bool getResidual(VarVectorCollection& residualv);
394 virtual void increaseSubspaceLimit(int newSize);
395
396 ergo_real getFreq(int i) const { return ritzVals[i]; }
397 void computeMoments(OneElOperator& dipx,
398 OneElOperator& dipy,
399 OneElOperator& dipz);
400 ergo_real getTransitionMoment2(int i) const { return transMoms2[i]; }
401};
402
403} /* End of namespace LR */
404#endif /* !defined(SLR_HEADER) */
E2Evaluator interface provides a way to perform a linear transformation of supplied transition densit...
Definition slr.h:181
virtual bool transform(const ergo_real *dmat, ergo_real *fmat)=0
virtual ~E2Evaluator()
Definition slr.h:184
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:1192
void computeMoments(OneElOperator &dipx, OneElOperator &dipy, OneElOperator &dipz)
Definition slr.cc:1226
ergo_real getTransitionMoment2(int i) const
Definition slr.h:400
ergo_real * ritzVals
recent ritz values in the subspace.
Definition slr.h:368
virtual bool getResidual(VarVectorCollection &residualv)
get residual of the eigenvalue problem.
Definition slr.cc:1058
ergo_real * last_ev
most recent eigenvectors in the reduced space
Definition slr.h:372
virtual void increaseSubspaceLimit(int newSize)
expands above the default limit
Definition slr.cc:1173
ergo_real getFreq(int i) const
Definition slr.h:396
int nConverged
number of already converged eigenvalues
Definition slr.h:371
virtual ergo_real getPreconditionerShift(int i) const
returns the preconditioning shift.
Definition slr.h:388
ergo_real * transMoms2
most recent SQUARED transition moments.
Definition slr.h:369
EigenSolver(int nbast, int nocc, const ergo_real *fock_matrix, const ergo_real *overlap, int n)
Definition slr.h:374
int nStates
number of excited states to compute
Definition slr.h:370
virtual ~EigenSolver()
Definition slr.h:383
VarVectorCollection Avects
vects and Avects members store the trial vectors and their transformed versions.
Definition slr.h:325
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:389
VarVector e2diag
approximation to the diagonal of E2 operator
Definition slr.h:286
static const int MVEC
default limit for subspace size
Definition slr.h:285
ergo_real * cmo
the MO coefficients.
Definition slr.h:328
VarVectorCollection vects
base vectors
Definition slr.h:317
virtual int getInitialGuess(VarVectorCollection &vecs)=0
Computes the initial vector the subspace is to be seeded with.
bool lintrans(E2Evaluator &e2, const VarVector &v, VarVector &Av) const
performs the linear transformation of the vector with E[2] operator.
Definition slr.cc:707
void getAvMinusFreqSv(ergo_real f, ergo_real *weights, VarVector &r)
Computes a vector built of base vectors with specified vectors.
Definition slr.cc:852
virtual bool getResidual(VarVectorCollection &residualv)=0
Computes the residual vector.
int nbast
number of basis functions
Definition slr.h:315
ergo_real * fdiag
the eigenvalues of the Fock matrix.
Definition slr.h:326
SmallMatrix eSub
E[2] matrix projected onto subspace.
Definition slr.h:288
void operToVec(OneElOperator &oper, VarVector &res) const
Transform square operator to the vector form.
Definition slr.cc:902
void mo2ao(int nbast, const ergo_real *mo, ergo_real *ao) const
Definition slr.cc:604
virtual void increaseSubspaceLimit(int newSize)
expands above the default limit
Definition slr.cc:735
int maxSubspaceSize
current subspace size limit.
Definition slr.h:283
virtual ~LRSolver()
Definition slr.h:253
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:804
ergo_real convThreshold
iterative method convergence threshold
Definition slr.h:282
bool solve(E2Evaluator &e, bool diskMode=false)
Solves the problem defined by the subclass.
Definition slr.cc:753
int nocc
number of occupied orbitals
Definition slr.h:316
void buildVector(const ergo_real *w, VarVector &full)
Build full fector from the reduced form.
Definition slr.cc:889
int subspaceSize
current subspace size
Definition slr.h:287
void projectOnSubspace(const VarVector &full, ergo_real *w)
Projects vector.
Definition slr.cc:876
ergo_real * xSub
solution vector projected onto subspace
Definition slr.h:290
virtual void addToSpace(VarVectorCollection &vecs, E2Evaluator &e2)
extends the subspace with v and its transformed vector Av.
Definition slr.cc:641
virtual ergo_real getPreconditionerShift(int i) const =0
returns the preconditioning shift.
void ao2mo(int nbast, const ergo_real *ao, ergo_real *mo) const
computes mo := cmo'*ao*cmo
Definition slr.cc:615
void computeExactE2Diag(E2Evaluator &e2)
Definition slr.cc:788
SmallMatrix sSub
S[2] matrix projected onto subspace.
Definition slr.h:289
void load_F_MO(ergo_real *fmat) const
Definition slr.cc:625
Abstract interface to a one electron operator.
Definition slr.h:216
virtual ~OneElOperator()
Definition slr.h:219
virtual void getOper(ergo_real *result)=0
VarVector rhs
RHS of the SOE.
Definition slr.h:339
virtual bool getResidual(VarVectorCollection &residualv)
get the residual of the set of linear equations.
Definition slr.cc:967
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:343
virtual void increaseSubspaceLimit(int newSize)
expands above the default limit
Definition slr.cc:996
virtual int getInitialGuess(VarVectorCollection &vecs)
returns the initial guess for the linear set of equations.
Definition slr.cc:926
ergo_real frequency
frequency for which the SOE is to be solved.
Definition slr.h:338
virtual void addToSpace(VarVectorCollection &vecs, E2Evaluator &e2)
extends the subspace with v and its transformed vector Av.
Definition slr.cc:1022
void setRHS(OneElOperator &op)
initializes the rhs field
Definition slr.cc:1008
virtual ~SetOfEqSolver()
Definition slr.h:350
virtual ergo_real getPreconditionerShift(int) const
returns the preconditioning shift.
Definition slr.h:351
ergo_real getPolarisability(OneElOperator &oper)
computes polarizability by contracting the response vector with specified operator
Definition slr.cc:1039
ergo_real * rhsSub
RHS vector projected onto subspace.
Definition slr.h:358
ergo_real multiplyXtimesVec(const VarVector &rhs)
multiplies current solution by some vector.
Definition slr.cc:950
ergo_real xTimesRHS
Definition slr.h:361
a class implementing dynamic resized two dimensional arrays.
Definition slr.h:223
SmallMatrix(int sz)
Definition slr.h:235
int nsize
Definition slr.h:225
~SmallMatrix()
Definition slr.h:236
void expand(int newSize)
increase the dimension of the matrix without losing the data.
Definition slr.cc:374
const RowProxy operator[](int i)
Definition slr.h:237
a collection of vectors, usually handled at once.
Definition slr.h:189
int nAllocated
Definition slr.h:194
void release()
Make sure there is space for one vector.
Definition slr.cc:330
static const char * tmpdir
Definition slr.h:211
~VarVectorCollection()
Definition slr.cc:302
void setDiskMode(bool x)
Definition slr.h:206
unsigned * ages
Definition slr.h:191
VarVector & operator[](int i)
Definition slr.cc:319
unsigned currentAge
Definition slr.h:192
void setSize(int sz)
Definition slr.cc:309
VarVector * vecs
Definition slr.h:190
int size() const
Definition slr.h:204
bool diskMode
Definition slr.h:195
VarVectorCollection(int nSize=0)
Definition slr.h:197
int nVecs
Definition slr.h:193
void releaseAll()
Release all vectors from the memory, saving if necessary.
Definition slr.cc:356
bool getDiskMode() const
Definition slr.h:205
template based proxy object that uses bool-valued policies to perform the assignments.
Definition slr.h:51
ergo_real scalar
Definition slr.h:54
VarVectorProxyOp(const VarVector &a, ergo_real s=1.0)
Definition slr.h:55
const VarVector & vec
Definition slr.h:53
Vector of variables parametrising the solution to the linear response equations.
Definition slr.h:62
ergo_real * v
vector coefficients
Definition slr.h:63
ergo_real & operator[](int i)
Definition slr.h:114
VarVector()
Definition slr.h:79
void setFromFull(int nbast, int nocc, const ergo_real *full_mat)
Definition slr.cc:404
VarVector & operator=(const VarVector &b)
Definition slr.h:120
int nvar
nvar := nocc*nvirt.
Definition slr.h:67
VarVector & operator=(ergo_real scalar)
Definition slr.h:115
void symorth(void)
Uses symmetric orthogonalization to orthonormalize itself (x y) with the swapped vector (y x).
Definition slr.cc:466
void setFull(int nbast, int nocc, ergo_real *full_mat) const
Definition slr.cc:433
void setSize(int n)
Definition slr.h:97
ergo_real * y() const
returns the Y part.
Definition slr.h:95
const ergo_real & operator[](int i) const
Definition slr.h:113
void print(const char *comment=NULL)
Definition slr.h:107
unsigned inMemory
valid representation in memory
Definition slr.h:69
VarVector & operator=(const VarVectorProxyOp< false, true > &proxy)
Definition slr.h:148
VarVector & operator=(const VarVectorProxyOp< false, false > &proxy)
Definition slr.h:133
VarVector(const VarVector &a)
Definition slr.h:71
char * fName
Present in secondary storage (i.e.
Definition slr.h:65
VarVector(int nbast, int nocc, const ergo_real *full_mat)
Creates a vector from a full matrix.
Definition slr.h:82
~VarVector()
Definition slr.h:87
ergo_real * x() const
return the X part of the vector.
Definition slr.h:94
int size() const
Definition slr.h:106
friend ergo_real operator*(const VarVector &a, const VarVector &b)
Definition slr.cc:101
unsigned onDisk
valid representation on disk
Definition slr.h:68
void save(const char *tmpdir)
Save the object.
Definition slr.cc:250
void load(const char *tmpdir)
Load the object to memory.
Definition slr.cc:219
void release(const char *tmpdir)
Releases the memory, saving if necessary.
Definition slr.cc:286
Definition slr.cc:90
Definition allocate.cc:39
Definition of the main floating-point datatype used; the ergo_real type.
double ergo_real
Definition realtype.h:69
Definition slr.h:227
ergo_real & operator[](int j) const
Definition slr.h:230
RowProxy(ergo_real *r)
Definition slr.h:229
ergo_real * toprow
Definition slr.h:228