-
Notifications
You must be signed in to change notification settings - Fork 2
/
crystal.py
121 lines (97 loc) · 4.35 KB
/
crystal.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
# Class for setting and getting crystal parameters, lattice vectors, and structure functions
from .materials import Material
from .constants import *
from .fmath import *
"""
Crystal class and methods
Inherits attributes from Material class
lattice_const: lattice constant in Angstroms
primitives: an array of the N primitive basis vectors in format [alpha0, alpha1,...]
for each alpha = [#,#,#] as a 3-list
a1, a2, a3: Bravais lattice vectors as 3-lists [#,#,#]
"""
class Crystal(Material):
def __init__(self, material, primitives, a1, a2, a3, volume=1.0, density=1.0,
fiducial_mass=1.0):
super().__init__(material_name=material, fiducial_mass=fiducial_mass, volume=volume, density=density)
self.a = self.lattice_const
self.alpha = self.a * np.array(primitives)
self.a0 = self.a * np.array([0,0,0])
self.a1 = self.a * np.array(a1)
self.a2 = self.a * np.array(a2)
self.a3 = self.a * np.array(a3)
self.basis = np.array([self.a0, self.a1, self.a2, self.a3])
self.b1 = 2*pi*cross(self.a2, self.a3) / (dot(self.a1, cross(self.a2, self.a3)))
self.b2 = 2*pi*cross(self.a3, self.a1) / (dot(self.a1, cross(self.a2, self.a3)))
self.b3 = 2*pi*cross(self.a1, self.a2) / (dot(self.a1, cross(self.a2, self.a3)))
self.norm = [0.0, 0.0, 1.0]
self.lx = power(volume, 1/3)
self.ly = power(volume, 1/3)
self.lz = power(volume, 1/3)
def r(self, n1, n2, n3):
return n1 * self.a1 + n2 * self.a2 + n3 * self.a3
def G(self, h, k, l):
return h * self.b1 + k * self.b2 + l * self.b3
def wavelength(self, h, k, l):
return 2*pi/sqrt(dot(self.G(h, k, l), self.G(h, k, l)))
def energy(self, h, k, l):
return HBARC * sqrt(dot(self.G(h, k, l), self.G(h, k, l)))
def miller(self, h, k, l):
return np.array([h, k, l])
def sfunc(self, h, k, l):
return (1+exp(-1j * dot(self.alpha[1], self.G(h, k, l)))) \
* np.sum([exp(-1j*dot(self.G(h, k, l), avec)) for avec in self.basis])
def SF2(self, h, k, l):
return np.real(self.sfunc(h,k,l) * np.conjugate(self.sfunc(h,k,l)))
# Namelist
cryslist = ["Ge", "Si", "NaI", "CsI", "TeO2"]
def get_crystal(name, volume):
if name not in cryslist:
print("Specified material not in library. Supported crystals:\n", cryslist)
return
if name == "Ge":
# Diamond cubic, germanium
GeAlpha0 = [0.0, 0.0, 0.0]
GeAlpha1 = [0.25, 0.25, 0.25]
Primitives = [GeAlpha0, GeAlpha1]
A1 = [0.0, 0.5, 0.5]
A2 = [0.5, 0.0, 0.5]
A3 = [0.5, 0.5, 0.0]
return Crystal(name, Primitives, A1, A2, A3, volume=volume)
if name == "Si":
# Diamond cubic, sodium iodide
Alpha0 = [0.0, 0.0, 0.0]
Alpha1 = [0.25, 0.25, 0.25]
Primitives = [Alpha0, Alpha1]
A1 = [0.0, 0.5, 0.5]
A2 = [0.5, 0.0, 0.5]
A3 = [0.5, 0.5, 0.0]
return Crystal(name, Primitives, A1, A2, A3, volume=volume)
if name == "NaI":
# Diamond cubic, sodium iodide
Alpha0 = [0.0, 0.0, 0.0]
Alpha1 = [0.5, 0.5, 0.5]
Primitives = [Alpha0, Alpha1]
A1 = [0.0, 0.5, 0.5]
A2 = [0.5, 0.0, 0.5]
A3 = [0.5, 0.5, 0.0]
return Crystal(name, Primitives, A1, A2, A3, volume=volume)
if name == "CsI":
# Diamond cubic, cesium iodide
Alpha0 = [0.0, 0.0, 0.0]
Alpha1 = [0.5, 0.5, 0.5]
Primitives = [Alpha0, Alpha1]
A1 = [0.0, 0.5, 0.5]
A2 = [0.5, 0.0, 0.5]
A3 = [0.5, 0.5, 0.0]
return Crystal(name, Primitives, A1, A2, A3, volume=volume)
if name == "TeO2":
# D4(422) space group
Primatives_O = [[0.998, 0.916, 1.480], [1.481, 3.396, 3.386], [1.399, 3.314, 0.426], [0.916, 0.998, 6.145],
[3.797, 3.879, 5.293], [3.314, 1.399, 7.199], [3.396, 1.481, 4.239], [3.879, 3.797, 2.332]]
Primitives_Te = [[0.093, 0.093, 0.000],[4.702, 4.702, 3.813], [2.304, 2.491 , 1.906], [2.491, 2.304, 5.719]]
A1 = [0.0, 0.5, 0.5]
A2 = [0.5, 0.0, 0.5]
A3 = [0.5, 0.5, 0.0]
return Crystal(name, Primitives_Te+Primatives_O, A1, A2, A3, volume=volume)
return