ergo
Main Page
Namespaces
Classes
Files
File List
File Members
slr.h
Go to the documentation of this file.
1
/* Ergo, version 3.2, a program for linear scaling electronic structure
2
* calculations.
3
* Copyright (C) 2012 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>
41
class
VarVectorProxyOp
{
42
public
:
43
const
VarVector
&
vec
;
44
ergo_real
scalar
;
45
explicit
VarVectorProxyOp
(
const
VarVector
& a,
ergo_real
s=1.0) :
vec
(a),
scalar
(s) {}
46
};
47
48
52
class
VarVector
{
53
ergo_real
*
v
;
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
77
~VarVector
() {
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]; }
103
VarVector
&
operator=
(
ergo_real
scalar) {
104
for
(
int
i=0; i<2*
nvar
; i++)
v
[i] = scalar;
105
onDisk
= 0;
106
return
*
this
;
107
};
108
VarVector
&
operator=
(
const
VarVector
& b) {
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
&
121
operator=
(
const
VarVectorProxyOp<false, false >
& proxy) {
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
&
136
operator=
(
const
VarVectorProxyOp<false, true >
& proxy) {
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,
161
const
VarVectorProxyOp<false,false>
& b);
162
friend
ergo_real
operator*
(
const
VarVector
&a,
163
const
VarVectorProxyOp<true,false>
& b);
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
/* ################################################################### */
177
class
VarVectorCollection
{
178
VarVector
*
vecs
;
179
unsigned
*
ages
;
180
unsigned
currentAge
;
181
int
nVecs
;
182
int
nAllocated
;
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
}
189
~VarVectorCollection
();
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
/* ################################################################### */
204
class
OneElOperator
{
205
public
:
206
virtual
void
getOper
(
ergo_real
*result) = 0;
207
virtual
~OneElOperator
() {}
208
};
209
211
class
SmallMatrix
{
212
ergo_real
*
mat
;
213
int
nsize
;
214
protected
:
215
struct
RowProxy
{
216
ergo_real
*
toprow
;
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
);
269
void
computeExactE2Diag
(
E2Evaluator
& e2);
270
ergo_real
convThreshold
;
271
int
maxSubspaceSize
;
272
protected
:
273
static
const
int
MVEC
= 200;
274
VarVector
e2diag
;
275
int
subspaceSize
;
276
SmallMatrix
eSub
;
277
SmallMatrix
sSub
;
278
ergo_real
*
xSub
;
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
;
305
VarVectorCollection
vects
;
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
:
313
VarVectorCollection
Avects
;
314
ergo_real
*
fdiag
;
316
ergo_real
*
cmo
;
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
{
326
ergo_real
frequency
;
327
VarVector
rhs
;
328
public
:
331
SetOfEqSolver
(
int
nbast
,
int
nocc
,
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
:
346
ergo_real
*
rhsSub
;
347
virtual
void
addToSpace
(
VarVectorCollection
& vecs,
E2Evaluator
& e2);
348
ergo_real
multiplyXtimesVec
(
const
VarVector
&
rhs
)
/* const */
;
349
ergo_real
xTimesRHS
;
350
};
351
352
353
/* ################################################################### */
355
class
EigenSolver
:
public
LRSolver
{
356
ergo_real
*
ritzVals
;
357
ergo_real
*
transMoms2
;
358
int
nStates
;
359
int
nConverged
;
360
ergo_real
*
last_ev
;
361
public
:
362
EigenSolver
(
int
nbast
,
int
nocc
,
363
const
ergo_real
*fock_matrix,
364
const
ergo_real
*overlap,
int
n)
365
:
LRSolver
(nbast, nocc, NULL, NULL),
366
ritzVals
(new
ergo_real
[
MVEC
]),
transMoms2
(new
ergo_real
[
MVEC
]),
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) */
source
slr.h
Generated on Fri Dec 21 2012 12:20:28 for ergo by
1.8.1.1