-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathGenerateStartingFunctions.py
175 lines (129 loc) · 4.74 KB
/
GenerateStartingFunctions.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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
import numpy as np
import numpy.linalg as LA
import math
# This file generates starting wavefunctions for a particular configuration.
# 1. Write functions analytically
# 2. Load lookup table from txt
# Hydrogen H2 Molecule from atomic wavefunctions
# Atomic units
m_e = 1.0 #* 9.11e-31 # electron mass
q_e = 1.0 #* 1.602e-19 # electron charge
hbar = 1.0 #* 1.0546e-34 # reduced planck constant
k_e = 1.0 #/ (4*np.pi * 8.854e-12) # Coulomb's constant (1/4\pi\epsilon_0)
a_B = 1.0 #* hbar*hbar / (k_e*m_e*q_e**2) # Bohr radius
# Define atom positions
ion_positions = np.zeros((2,3))
ion_sep = 1.0
R1 = np.array([-0.5*ion_sep,0,0])*a_B
R2 = np.array([0.5*ion_sep,0,0])*a_B
ion_positions[0,:] = R1
ion_positions[1,:] = R2
ion_charges = np.array([1.0,1.0])
# Spherical Harmonics, indexed by l, m
Y00 = (4.0*math.pi)**(-0.5)
def Y1all(rvec, Z=1.0):
rmag = np.sqrt(np.sum(rvec*rvec,1))
return np.sqrt(0.75/np.pi) * rvec / np.array([rmag,rmag,rmag]).T
def Y1x(rvec, Z=1.0):
rmag = np.sqrt(np.sum(rvec*rvec,1))
return np.sqrt(0.75/np.pi) * rvec[:,0] / rmag
def Y1y(rvec, Z=1.0):
rmag = np.sqrt(np.sum(rvec*rvec,1))
return np.sqrt(0.75/np.pi) * rvec[:,1] / rmag
def Y1z(rvec, Z=1.0):
rmag = np.sqrt(np.sum(rvec*rvec,1))
return np.sqrt(0.75/np.pi) * rvec[:,2] / rmag
# Radial functions
# Note: this function requires that position coordinates are along dimension 1 of the array (not 0)
def R10(rvec, Z=1.0):
r = np.sqrt(np.sum(rvec*rvec,1))
a = a_B/Z
return 2*a**(-1.5)*np.exp(-r/a)
def R20(rvec, Z=1.0):
a = 2*a_B/Z
rb = np.sqrt(np.sum(rvec*rvec,1)) / a
return 2*a**(-1.5)*(1.0-rb)*np.exp(-rb)
def R21(rvec, Z=1.0):
a = 2*a_B/Z
rb = np.sqrt(np.sum(rvec*rvec,1)) / a
return 2*3**(-.5) * a**(-3.5) * rb * np.exp(-rb)
def Laplacian_R10(rvec, Z=1.0):
a = a_B/Z
rb = np.sqrt(np.sum(rvec*rvec,1)) / a
#return 2*a**(-3.5)*np.exp(-r/a)
return 2*a**(-3.5)*(1.0 - 2.0/rb) * np.exp(-rb)
def Laplacian_R20(rvec, Z=1.0):
a = 2*a_B/Z
rb = np.sqrt(np.sum(rvec*rvec,1)) / a # change of vars makes formula simpler: rb = r/a = r*Z / 2*a_B
#return 4*(2*a_B)**(-3.5)*(1.5-(r/(4.0*a_B)))*np.exp(-2.0*r/a_B)
#return (1.0/(2.0*math.sqrt(2.0)))*a_B**(-3.5)*(1.5-(r/(4.0*a_B)))*np.exp(-2.0*r/a_B)
return -2*a**(-3.5)/rb * (rb-4)*(rb-1)*np.exp(-rb)
####################################
#
# S orbitals
def psi_1s(e_pos_vec,i_pos,Z=1.0):
return Y00 * R10(e_pos_vec - i_pos, Z)
#def psi_s1(rvec):
# return Y00 * R10(rvec-R1)
#
#def psi_s2(rvec):
# return Y00 * R10(rvec-R2)
def psi_2s(e_pos_vec, i_pos, Z=1.0):
return Y00 * R20(e_pos_vec - i_pos, Z)
def psi_2p_all(e_pos_vec, i_pos, Z=1.0):
return Y1all(e_pos_vec - i_pos) * R21(e_pos_vec - i_pos, Z)
def psi_2px(e_pos_vec, i_pos, Z=1.0):
return Y1x(e_pos_vec - i_pos)[0] * R21(e_pos_vec - i_pos, Z)
def psi_2py(e_pos_vec, i_pos, Z=1.0):
return Y1y(e_pos_vec - i_pos)[1] * R21(e_pos_vec - i_pos, Z)
def psi_2pz(e_pos_vec, i_pos, Z=1.0):
return Y1z(e_pos_vec - i_pos)[2] * R21(e_pos_vec - i_pos, Z)
# Laplacian of S orbitals
def Lpsi_1s(e_pos_vec, i_pos):
return Y00 * Laplacian_R10(e_pos_vec - i_pos) * (-hbar*hbar*0.5/m_e)
def Lpsi_s1(rvec):
return Y00 * Laplacian_R10(rvec-R1) # TODO
def Lpsi_s2(rvec):
return Y00 * Laplacian_R10(rvec-R2) # TODO
# Retrieve the functions
def getH2Functions():
return np.array([psi_s1, psi_s2])
# Retrieve the Laplacians of the wavefunctions (can be defined analytically)
def getH2Laplacians():
return np.array([Lpsi_s1, Lpsi_s2])
class H_atom:
Z = 1.0
i_pos = np.zeros(3)
def __init__(self,pos=np.array([0,0,0])):
self.i_pos = pos
def setPosition(self, pos):
self.i_pos = pos
def psi_1s(self, e_pos_vec):
return Y00 * R10(e_pos_vec - self.i_pos, self.Z)
class Atom:
Z = 1.0
i_pos = np.zeros(3)
last_e_vec = np.zeros(3)
last_2p_vec = np.zeros(3)
def __init__(self, pos=np.array([0,0,0]), Z=1.0):
self.i_pos = pos
self.Z = float(Z)
def setPosition(self, pos):
self.i_pos = pos
def psi_1s(self, e_pos_vec):
return Y00 * R10(e_pos_vec - self.i_pos, self.Z)
def psi_2s(self, e_pos_vec):
return Y00 * R20(e_pos_vec - self.i_pos, self.Z)
def psi_2p_all(self, e_pos_vec):
return Y1all(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos, self.Z)
def psi_2px(self, e_pos_vec):
#if (self.last_e_vec == e_pos_vec).all():
# return self.last_2p_vec[0]
#else:
# self.last_e_vec = e_pos_vec.copy()
# self.last_2p_vec = Y1all(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos, self.Z)
return Y1x(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos,self.Z)
def psi_2py(self, e_pos_vec):
return Y1y(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos,self.Z)
def psi_2pz(self, e_pos_vec):
return Y1z(e_pos_vec - self.i_pos) * R21(e_pos_vec - self.i_pos,self.Z)