-
Notifications
You must be signed in to change notification settings - Fork 1
/
ex3.py
143 lines (117 loc) · 3.18 KB
/
ex3.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
from numpy import *
import copy
def energy(spatial_index):
if spatial_index == 1:
return 1.0
elif ((spatial_index == 2) | (spatial_index == 3)):
return 2.0
elif ((spatial_index == 4) | (spatial_index == 5) | (spatial_index == 6)):
return 3.0
def delta(i,j):
if(i == j):
return 1
else:
return 0
def density_matrix(U,L,N):
D = zeros((L/2,L/2))
for r in range(0,L/2):
for s in range(0,L/2):
D[s,r] = 0.0
for j in range(0,N/2):
D[s,r] += 2.0*U[s,j]*conjugate(U[r,j])
return D
def QRPS(q,r,p,s,w):
key1 = str(q+1)+str(r+1) + str(p+1) + str(s+1)
key2 = str(q+1)+str(r+1) + str(s+1) + str(p+1)
if((w.has_key(key1) == True) and (w.has_key(key2) == True)):
qrps = w[key1]-0.5*w[key2]
elif ((w.has_key(key1) == True) and (w.has_key(key2) == False)):
qrps = w[key1]
elif ((w.has_key(key1) == False) and (w.has_key(key2) == True)):
qrps = -0.5*w[key2]
else:
qrps = 0.0
return qrps
def compute_FockMatrix(D,L,w):
FD = zeros((L/2,L/2))
for q in range(0,L/2):
for p in range(0,L/2):
FD[q,p] = delta(q+1,p+1)*energy(p+1)
for r in range(0,L/2):
for s in range(0,L/2):
FD[q,p] += D[s,r]*QRPS(q,r,p,s,w)
return FD
def diag_FD(FD):
eps, U = linalg.eig(FD)
eps, U = bubbleSort(eps,U)
return eps, U
def bubbleSort(alist,eig_vec):
sorted_alist = copy.copy(alist)
sorted_eigvec = copy.copy(eig_vec)
for passnum in range(len(sorted_alist)-1,0,-1):
for i in range(passnum):
if sorted_alist[i]>sorted_alist[i+1]:
temp = sorted_alist[i]
temp2 = copy.copy(sorted_eigvec[:,i])
sorted_alist[i] = sorted_alist[i+1]
sorted_alist[i+1] = temp
sorted_eigvec[:,i] = copy.copy(sorted_eigvec[:,i+1])
sorted_eigvec[:,i+1] = temp2
return sorted_alist, sorted_eigvec
def check_F_Hermitian(F,L):
isHermitian = True
for p in range(0,L/2):
for q in range(0,L/2):
if (F[p,q] != conjugate(F[q,p])):
isHermitian = False
break
return isHermitian
def HF_Energy(eps, U, D, w, N, L):
ERHF = 0.0
for i in range(0,N/2):
ERHF += eps[i]
ERHF = 2.0*ERHF
for i in range(0,N/2):
for p in range(0,L/2):
for q in range(0,L/2):
tmp = 0.0
for s in range(0,L/2):
for r in range(0,L/2):
tmp += D[s,r]*QRPS(q,r,p,s,w)
ERHF -= conjugate(U[q,i])*tmp*U[p,i]
return ERHF
def computeHartreeFockSolution(L,N,w):
epsilon = 1*10**(-16)
U = identity(L/2)
D = density_matrix(U,L,N)
FD = compute_FockMatrix(D,L,w)
eps_old, U = diag_FD(FD)
max_iters = 30
ERHF = HF_Energy(eps_old,U,D,w,N,L)
print "ERHF = %.16f, from initial guess." % ERHF
for k in range(1,max_iters):
D = density_matrix(U,L,N)
FD = compute_FockMatrix(D,L,w)
eps_new, U = diag_FD(FD)
ERHF = HF_Energy(eps_new,U,D,w,N,L)
#print eps_new
print "ERHF = %.18f, |eps_new - eps_old| = %g" % (ERHF,max(abs(eps_new-eps_old)))
if(max(abs(eps_new-eps_old)) < epsilon):
eps_old = eps_new
break
eps_old = eps_new
return FD, U, eps_old, ERHF
inFile = open('coulomb.dat','r')
w = {}
for line in inFile:
tmp = line.split()
key = tmp[0] + tmp[1] + tmp[2] + tmp[3]
val = float(tmp[4])
w[key] = val
L = 12
N = 2
F, U, eps_old, ERHF = computeHartreeFockSolution(L,N,w)
for i in range(0,N/2):
for a in range(N/2,L/2):
D_ia = F[i,i]-F[a,a]
print D_ia