forked from Masterluke87/psixas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
kshelper.py
145 lines (122 loc) · 3.97 KB
/
kshelper.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
# -*- coding: utf-8 -*-
"""
Created on Sun Aug 19 18:30:11 2018
@author: luke
"""
import numpy as np
import time
def diag_H(H, A):
Hp = A.dot(H).dot(A)
e, C2 = np.linalg.eigh(Hp)
C = A.dot(C2)
return (C,e)
class Timer(object):
def __init__(self):
self.entries = {}
def addStart(self,entry):
if (entry in self.entries):
self.entries[entry]["Start"] = time.time()
else:
self.entries[entry] = { }
self.entries[entry]["Start"] = time.time()
def addEnd(self,entry):
if (entry in self.entries):
self.entries[entry]["End"] = time.time()
else:
self.entries[entry] = { }
self.entries[entry]["End"] = time.time()
def getTime(self,entry):
return (self.entries[entry]["End"] - self.entries[entry]["Start"])
def printAlltoFile(self,filename):
S = "Timings:\n"
for x in self.entries:
S += "{:^10} | {:5.2f}s \n".format(x,self.getTime(x))
S+="\n"
open(filename,"a").write(S)
class DIIS_helper(object):
"""
A helper class to compute DIIS extrapolations.
Notes
-----
Equations taken from [Sherrill:1998], [Pulay:1980:393], & [Pulay:1969:197]
Algorithms adapted from [Sherrill:1998] & [Pulay:1980:393]
"""
def __init__(self, max_vec=6):
"""
Intializes the DIIS class.
Parameters
----------
max_vec : int (default, 6)
The maximum number of vectors to use. The oldest vector will be deleted.
"""
self.error = []
self.vector = []
self.max_vec = max_vec
def add(self, state, error):
"""
Adds a set of error and state vectors to the DIIS object.
Parameters
----------
state : array_like
The state vector to add to the DIIS object.
error : array_like
The error vector to add to the DIIS object.
Returns
------
None
"""
error = np.array(error)
state = np.array(state)
if len(self.error) > 1:
if self.error[-1].shape[0] != error.size:
raise Exception("Error vector size does not match previous vector.")
if self.vector[-1].shape != state.shape:
raise Exception("Vector shape does not match previous vector.")
self.error.append(error.ravel().copy())
self.vector.append(state.copy())
def extrapolate(self):
"""
Performs the DIIS extrapolation for the objects state and error vectors.
Parameters
----------
None
Returns
------
ret : ndarray
The extrapolated next state vector
"""
# Limit size of DIIS vector
diis_count = len(self.vector)
if diis_count == 0:
raise Exception("DIIS: No previous vectors.")
if diis_count == 1:
return self.vector[0]
if diis_count > self.max_vec:
# Remove oldest vector
del self.vector[0]
del self.error[0]
diis_count -= 1
# Build error matrix B
B = np.empty((diis_count + 1, diis_count + 1))
B[-1, :] = -1
B[:, -1] = -1
B[-1, -1] = 0
for num1, e1 in enumerate(self.error):
B[num1, num1] = np.vdot(e1, e1)
for num2, e2 in enumerate(self.error):
if num2 >= num1: continue
val = np.vdot(e1, e2)
B[num1, num2] = B[num2, num1] = val
# normalize
B[abs(B) < 1.e-14] = 1.e-14
B[:-1, :-1] /= np.abs(B[:-1, :-1]).max()
# Build residual vector
resid = np.zeros(diis_count + 1)
resid[-1] = -1
# Solve pulay equations
ci = np.dot(np.linalg.pinv(B), resid)
# combination of previous fock matrices
V = np.zeros_like(self.vector[-1])
for num, c in enumerate(ci[:-1]):
V += c * self.vector[num]
return V