forked from erwinjr/erwinjr
-
Notifications
You must be signed in to change notification settings - Fork 1
/
QCLayers.py
1230 lines (1132 loc) · 55.7 KB
/
QCLayers.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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python2
# -*- coding:utf-8 -*-
# ===========================================================================
# ErwinJr is a simulation program for quantum semiconductor lasers.
# Copyright (C) 2012 Kale J. Franz, PhD
# Copyright (C) 2017 Ming Lyu (CareF)
#
# A portion of this code is Copyright (c) 2011, California Institute of
# Technology ("Caltech"). U.S. Government sponsorship acknowledged.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# ===========================================================================
# TODO:
# material related codes should be moved to MaterialConstants
# try use dict type for substrate restriction on material
# Add solve well (helping tight binding model design
# Performance improve: lo_transition_rate and solve_psi:
# parallel support: Done
# reduce reductant calculation
# Graphic card support?
# Try matrix eigen solver in solve_psi - It's O(n^3), compared to what we
# have for O(n^2) (or O(mn) where m for E resolution and n for x)
# replace CLIB by Cython
from __future__ import division
__LOG__ = False
__DEBUG__ = 1
__USE_CLIB__ = True
__MORE_INTERPOLATION__ = True # One more time interpolation for eigen solver
__MULTI_PROCESSING__ = True
import copy
import sys
import numpy as np
from numpy import sqrt, exp, pi
from scipy import interpolate
from settings import (wf_scale, psi_scale, wf_min_height, pretty_plot_factor,
plot_decimate_factor, phonon_integral_factor)
import MaterialConstantsDict
cst = MaterialConstantsDict.MaterialConstantsDict()
if __LOG__:
import pickle
logcount = 0
if __DEBUG__ >= 1:
pass
# from time import time
# TODO: replace CLIB by Cython
if __USE_CLIB__:
# from ctypes import *
import ctypes as ct
if __MULTI_PROCESSING__:
if sys.platform in ('linux2', 'darwin', 'cygwin'):
cQ = np.ctypeslib.load_library('cQCLayersMP', '.')
elif sys.platform == 'win32':
cQ = ct.CDLL('cQCLayersMP.dll')
else:
if sys.platform in ('linux2', 'darwin', 'cygwin'):
cQ = np.ctypeslib.load_library('cQCLayers', '.')
elif sys.platform == 'win32':
cQ = ct.CDLL('cQCLayers.dll')
# ===========================================================================
# Global Variables
# ===========================================================================
from scipy.constants import (e as e0, epsilon_0 as eps0,
electron_mass as m0, c as c0)
from scipy.constants import h, hbar
ANG = 1e-10 # angstrom to meter
KVpCM = 1e5 # KV/cm to V/m
meV = 1e-3 # meV to eV
INV_INF = 1e-20 # for infinit small decay rate (ns-1)
PAD_HEAD = 100 # width padded in the head of the given region for basis solver
PAD_TAIL = 30
# ===========================================================================
# Reference
# [0]Kale Franz's thesis
# [1]Handbook of Optics, Vol.2, ISBN: 0070479747
# [2]Van de Walle C G. Band lineups and deformation potentials in the
# model-solid theory[J]. Physical review B, 1989, 39(3): 1871.
# [3]Peter Qiang Liu's thesis
# ===========================================================================
# for In0.53Ga0.47As, EcG = 0.22004154
# use this as a zero point baseline
bandBaseln = 0.22004154
class QCLayers(object):
"""Class for QCLayers
Member variables:
parameters for each layer, np.array type, with len = No. of layers:
layerWidth -int in pixel (xres per pixel), width of each layer
layerBarriers -boolean(TBD), if the layer is barrier or not
layerARs -boolean(TBD), if the layer is active region or not
only affect basis solver (negelet some coupling
layerMaterials -int(TBD), label of material, depending on
substrate, the material is defined in erwinjr.pyw
layerDopings -Doping per volumn in unit 1e17 cm-3
layerDividers -??? seems not used
xres: position resolution, in angstrom
vertRes: vertical/energy resolution, in meV
EField: external (static) electrical field, in kV/cm = 1e5 V/m
layerSelected: (for GUI) a label which layer is selected in GUI,
with default -1 indicating non selected.
repeats: (int) is the number of repeat times for the given structure
solver: ??? seems not used
Temperature: Temperature of the device, affecting material property
seems not used
TempFoM: ??? seems not used (Figure of Merit?
diffLength: ??? seems not used
substrate: The substrate material for the device, which determined
the well and barrier material
substrate | well | barrier
InP | In_xGa_{1-x}As | Al_{1-x}In_xAs
GaAs | Al_xGa_{1-x}As | Al_xGa_{1-x}As
GaSb | InAs_ySb_{1-y} | Al_xGa_{1-x}Sb
GaN (TBD)
basisARInjector & basisInjectorAR: where should the layer separate
for basis solver
moleFrac: mole fraction for each possible layer material, in format
[well, barrier]*4
"""
def __init__(self):
self.layerWidth = np.array([1, 1]) # pix
self.layerBarriers = np.array([0, 0]) # boolean
self.layerARs = np.array([0, 0]) # boolean
self.layerMaterials = np.array([1, 1]) # label
self.layerDopings = np.array([0., 0.]) # 1e17 cm-3
self.layerDividers = np.array([0, 0]) # Basis solver divider
self.xres = 0.5 # angstrom
self.EField = 0 # kV/cm = 1e5 V/m
self.layerSelected = -1 # int
self.vertRes = 0.5 # meV
self.repeats = 2 # repeats n times for the given structure
self.description = ""
self.solver = "SolverH" # ?
self.Temperature = 300
self.TempFoM = 300 # ?
self.diffLength = 0 # ?
self.basisARInjector = True
self.basisInjectorAR = True
# self.designByAngs = True
# self.designByML = False
self.substratesList = ('InP', 'GaAs', 'GaSb', 'GaN')
self.substrate = 'InP'
self.moleFrac = [0.53, 0.52, 0.53, 0.52, 0.53, 0.52, 0.53, 0.52]
self.update_alloys()
self.update_strain()
self.populate_x()
def set_xres(self, res):
for n in range(self.layerWidth.size):
self.layerWidth[n] = int(np.round(
self.layerWidth[n] * self.xres / res))
self.xres = res
def populate_x(self):
"""Extend layer information to position functions
Layer data: layerWidth
with len = # of layers and each value repr. a layer
Position data (OUTPUT/update member variables):
xPoints
position grid
xBarriers: from layerBarriers, is barrier layer
should be boolean (TBD)
xARs: from layerARs and xVc, xVc if is active region
(layerARs == 1) otherwise np.NaN
xMaterials: from layerMaterials label/index of material
should be int starting from 0 (TBD)
xDopings: from layerDopings, doping per volumn
xLayerNums
at xPoints[q] it's xLayerNums[q]-th layer
xLayerSelected: from layerSelected and xVc, xVc if it's
the layer indicated by layerSelected, otherwise NaN
"""
# print "-----debug----- QCLayers populate_x called"
# print self.layerBarriers
self.xPoints = self.xres * np.arange(0, self.layerWidth.sum())
layerNumCumSum = np.concatenate(([0], self.layerWidth.cumsum()))
self.xBarriers = np.zeros(self.xPoints.shape)
self.xARs = np.zeros(self.xPoints.shape)
self.xMaterials = np.zeros(self.xPoints.shape)
self.xDopings = np.zeros(self.xPoints.shape)
self.xLayerNums = np.zeros(self.xPoints.shape)
# extend layer data for all xpoints
for q in xrange(0, self.layerWidth.size):
self.xBarriers[layerNumCumSum[q]:
layerNumCumSum[q + 1]] = self.layerBarriers[q]
if self.layerARs[q] == 1:
self.xARs[layerNumCumSum[q] - 1:
layerNumCumSum[q + 1] + 1] = 1
self.xMaterials[layerNumCumSum[q]:
layerNumCumSum[q + 1]] = self.layerMaterials[q]
self.xDopings[layerNumCumSum[q]:
layerNumCumSum[q + 1]] = self.layerDopings[q]
self.xLayerNums[layerNumCumSum[q]:
layerNumCumSum[q + 1]] = q
# duplicate layer based on user input repeats
# repeats = self.repeats
if self.repeats >= 2:
self.xPoints = np.arange(
0, self.xres * (self.layerWidth.sum() +
self.layerWidth[1:].sum() * (
self.repeats - 1)),
self.xres)
self.xBarriers = np.concatenate((self.xBarriers, np.tile(
self.xBarriers[layerNumCumSum[1]:], self.repeats - 1)))
self.xARs = np.concatenate((self.xARs, np.tile(
self.xARs[layerNumCumSum[1]:], self.repeats - 1)))
self.xMaterials = np.concatenate((self.xMaterials, np.tile(
self.xMaterials[layerNumCumSum[1]:], self.repeats - 1)))
self.xDopings = np.concatenate((self.xDopings, np.tile(
self.xDopings[layerNumCumSum[1]:], self.repeats - 1)))
self.xLayerNums = np.concatenate((self.xLayerNums, np.tile(
self.xLayerNums[layerNumCumSum[1]:], self.repeats - 1)))
# this hack is needed because sometimes
# self.xPoints is one element too big
self.xPoints = self.xPoints[0:self.xBarriers.size]
self.update_strain()
# Following are equiv. elec potential for different bands
# external field is included
# xVX, xVL, xVLH and xVSO are used for checking if there's indrect
# bandgap, s.t. we can prevent its effect
self.xVc = np.zeros(self.xPoints.size)
self.xVX = np.zeros(self.xPoints.size)
self.xVL = np.zeros(self.xPoints.size)
self.xVLH = np.zeros(self.xPoints.size)
self.xVSO = np.zeros(self.xPoints.size)
for MLabel in range(1, 5):
indx = np.nonzero(self.xMaterials == MLabel)[0]
if indx.size != 0:
material = np.where(self.xBarriers[indx] == 1,
MLabel * 2 - 1, (MLabel - 1) * 2)
self.xVc[indx] = (self.EcG[material] - self.xPoints[indx] *
ANG * self.EField * KVpCM)
self.xVX[indx] = (self.EcX[material] - self.xPoints[indx] *
ANG * self.EField * KVpCM)
self.xVL[indx] = (self.EcL[material] - self.xPoints[indx] *
ANG * self.EField * KVpCM)
self.xVLH[indx] = (self.EvLH[material] - self.xPoints[indx] *
ANG * self.EField * KVpCM)
self.xVSO[indx] = (self.EvSO[material] - self.xPoints[indx] *
ANG * self.EField * KVpCM)
# make array to show selected layer in mainCanvas
try:
self.xLayerSelected = np.zeros(self.xPoints.shape) * np.NaN
layerSelected = self.layerSelected
if layerSelected != -1:
if layerSelected == 0:
# row for first layer is selected
for repeat in xrange(1, self.repeats + 1):
base = layerNumCumSum[-1] * (repeat - 1)
self.xLayerSelected[
base + layerNumCumSum[layerSelected]:
base + layerNumCumSum[layerSelected + 1] + 1
] = self.xVc[
base + layerNumCumSum[layerSelected]:
base + layerNumCumSum[layerSelected + 1] + 1
]
elif self.layerSelected == self.layerWidth.size:
# last (blank) layer row is selected
pass
else:
for repeat in xrange(1, self.repeats + 1):
base = np.sum(self.layerWidth[1:]) * (repeat - 1)
self.xLayerSelected[
base + layerNumCumSum[layerSelected] - 1:
base + layerNumCumSum[layerSelected + 1] + 1
] = self.xVc[
base + layerNumCumSum[layerSelected] - 1:
base + layerNumCumSum[layerSelected + 1] + 1
]
except IndexError:
# index error happens in SolveBasis when the selected layer is
# greater than the number of layers in the solve swath
# however, xLayerSelected is not used for the SolveBasis function
# print ("Index Error for layer selection at function"
# "qclayer.populate_x")
pass
self.xARs[np.nonzero(self.xARs == 0)[0]] = np.NaN
self.xARs *= self.xVc
def populate_x_band(self):
"""Extend layer information to position functions for band parameter
OUTPUT/update member variables):
xEg, xMc, xESO, xEp, xF,
whose value is determeined by the layer material
TODO: tell why we need populate_x and populate_x_band
"""
# print "------debug------- QCLayers populate_x_band called"
# Following parameters can be looked up in cQCLayers.c
self.xEg = np.zeros(self.xPoints.size)
self.xMc = np.zeros(self.xPoints.size) # Seems not to be used
self.xESO = np.zeros(self.xPoints.size)
self.xEp = np.zeros(self.xPoints.size)
self.xF = np.zeros(self.xPoints.size)
for MLabel in range(1, 5):
indx = np.nonzero(self.xMaterials == MLabel)[0]
if indx.size != 0:
material = np.where(self.xBarriers[indx] == 1,
MLabel * 2 - 1, (MLabel - 1) * 2)
self.xEg[indx] = self.EgLH[material]
self.xMc[indx] = self.me[material]
self.xESO[indx] = self.ESO[material]
self.xEp[indx] = self.Ep[material]
self.xF[indx] = self.F[material]
def update_alloys(self): # c is a Material_Constant class instance
""" update material parameter for the alloy used.
(Always followed by update_strain)
OUTPUT/update member variable:
all parameters listed in variables: np.array with len=numMaterials
labeled by sequence [well, barrier]*4
self.numMaterials: Number of differenc types of material
supported
self.epsrho: ???
"""
variables = (
'EgG', 'EgL', 'EgX', 'VBO', 'DSO', # unit eV
'me0', # seems not used
'acG', 'acL', 'acX', # Pikus-Bir interaction parameter
'Ep', 'F', # effective mass parameter, unit eV (Ep) and 1 (F)
'XiX', # strain correction to band at X point, unit eV
'b', 'av', 'alG', # strain correction to bands at Gamma, unit eV
'beG', 'alL', # Varsh correction
# 'beL', 'alX', 'beX', # seems not used
'epss', 'epsInf', # static and high-freq permitivity
'hwLO', # LO phonon energy, unit eV
'alc', # lattice const, unit angstrom
'c11', 'c12' # elestic stiffness constants
)
# print "----debug--- substrate is "+self.substrate
# substrate restriction on layer material,
# see doc string of QCLayers class
# Material are labeled by sequence [well, barrier]*4
if self.substrate == 'InP':
self.numMaterials = 8
self.Mat1 = ['InAs'] * 8
self.Mat2 = ['GaAs', 'AlAs'] * 4
MatCross = ['InGaAs', 'AlInAs'] * 4
elif self.substrate == 'GaAs':
self.numMaterials = 8
self.Mat1 = ['AlAs'] * 8
self.Mat2 = ['GaAs'] * 8
MatCross = ['AlGaAs'] * 8 # Note EgG_AlGaAs's moleFrac deps
elif self.substrate == 'GaSb':
self.numMaterials = 8
self.Mat1 = ['InAs', 'AlSb'] * 4
self.Mat2 = ['InSb', 'GaSb'] * 4
MatCross = ['InAsSb', 'AlGaSb'] * 4
# Note EgG's bowing moleFrac deps
else:
raise TypeError('substrate selection not allowed')
for item in variables:
setattr(self, item, np.empty(self.numMaterials))
para = getattr(self, item)
for n in range(self.numMaterials):
para[n] = getattr(cst[MatCross[n]],
item+'f')(self.moleFrac[n])
# para[n] = (self.moleFrac[n] *
# getattr(cst[self.Mat1[n]], item) +
# (1 - self.moleFrac[n]) *
# getattr(cst[self.Mat2[n]], item))
# if MatCross[n] in cst and hasattr(cst[MatCross[n]], item):
# # bowing parameter
# para[n] -= (self.moleFrac[n] * (1 - self.moleFrac[n]) *
# getattr(cst[MatCross[n]], item))
# See MaterialConstantsDict.py...
# if self.substrate == 'GaAs':
# for n in range(self.numMaterials):
# EgG_AlGaAs = cst['AlGaAs'].EgG + 1.310 * self.moleFrac[n]
# self.EgG[n] = (
# self.moleFrac[n] * cst['AlAs'].EgG +
# (1 - self.moleFrac[n]) * cst['GaAs'].EgG -
# self.moleFrac[n] * (1 - self.moleFrac[n]) * EgG_AlGaAs)
# elif self.substrate == 'GaSb':
# for n in range(1, self.numMaterials, 2):
# EgG_AlGaSb = cst['AlGaSb'].EgG + 1.22 * self.moleFrac[n]
# self.EgG[n] = (
# self.moleFrac[n] * cst['AlSb'].EgG +
# (1 - self.moleFrac[n]) * cst['GaSb'].EgG -
# self.moleFrac[n] * (1 - self.moleFrac[n]) * EgG_AlGaSb)
# set this once the others are set ???
self.epsrho = 1 / (1 / self.epsInf - 1 / self.epss)
def update_strain(self): # c is a Material_Constant class instance
"""Update strain and strain related parameters inside each layers
(Always called after update_alloys)
OUTPUT/update member variables:
(all below are np.array with len=numMaterials)
(Material are labeled by sequence [well, barrier]*4)
self.a_parallel: lattice const. within/parallel to the layer plane
self.eps_parallel: strain tensor within/parallel to the layer plane
self.a_perp: lattice const. perpendicular to the layer plane
self.eps_perp: strain tensor perpendicular to the layer plane
self.MaterialWidth: total width of a each material
self.netStrain: spacial average of eps_perp in unit of percentage
self.avghwLO: spacial average of hwLO in unit of eV
self.MLThickness: monolayer thickness? shown in GUI as
layerWidth/MLThickness. This is an average number of
layers and actually the edge is rough
self.Pec, self.Pe, self.Qe, self.Varsh: correction terms on bands,
See Kales's thesis, sec2
self.ESO: spin-orbit splitting, including strain correction
self.EgLH, self.EgSO: band bottom/top at Gamma Epoints respect to
conduction band
self.me: effective mass ignoring energy dependence, in unit m0
only used as a backup and for broadening_energy
may be deleted for further update
self.EcG, self.EcL, self.EcX: conduction band bottom at
Gamma, L and X points, respect to a give
baseline
self.EvLH, self.EvSO: valence band (LH/SO) top at Gamma point
(EcL, EcX, EvLH, EvSO are only used for plotting?)
"""
if self.substrate in cst.substrateSet:
self.a_parallel = cst[self.substrate].alc
# parallel littice constant depends on substrate
else:
raise TypeError('substrate selection not allowed')
# [2]Walle eqn 1b
self.eps_parallel = self.a_parallel / self.alc - 1
# [2]Walle eqn 2a and 4a
self.a_perp = self.alc * (1 - 2 * self.c12 / self.c11 *
self.eps_parallel)
# [2]Walle eqn 2b
self.eps_perp = self.a_perp / self.alc - 1
# = -2*self.c12/self.c11*self.eps_parallel
# total width of different material
self.MaterialWidth = np.zeros(self.numMaterials)
for i in range(4):
# Note that material are labeled by sequence [well, barrier]*4
# [1:] because first layer doesn't count
# indx = np.nonzero(self.layerMaterials[1:] == i+1)[0]
# # print i+1, self.layerMaterials
# (BUG FIXED: self.layerMaterials[1:] results in material index
# mismatch by 1)
# self.layerWidth includes an extra layer to promise first=last
indx = (self.layerMaterials == i + 1)
indx[0] = False # s.t. 1st layer doesn't count
self.MaterialWidth[2 * i + 1] = self.xres * np.sum(
self.layerWidth[indx] * self.layerBarriers[indx])
self.MaterialWidth[2 * i] = self.xres * np.sum(
self.layerWidth[indx]) - self.MaterialWidth[2 * i + 1]
# print "------debug-----", np.sum(self.MaterialWidth)
self.netStrain = 100 * np.sum(
self.MaterialWidth * self.eps_perp
) / np.sum(self.MaterialWidth) # in percentage
self.avghwLO = np.sum(self.MaterialWidth * self.hwLO
) / np.sum(self.MaterialWidth)
self.MLThickness = np.zeros(self.layerMaterials.size)
for n, (MLabel, BLabel) in enumerate(zip((1, 1, 2, 2, 3, 3, 4, 4),
(0, 1) * 4)):
# MLThickness is monolayer thickness of the material
self.MLThickness[
(self.layerMaterials == MLabel) &
(self.layerBarriers == BLabel)
] = self.a_perp[n] / 2.0
# Pikus-Bir interaction correction to bands offset,
# According to Kale's, Eq.(2.14),
# Pec for \delta E_{c} and Pe for \delta E_{v}
self.Pec = (2 * self.eps_parallel + self.eps_perp) * (self.acG)
self.Pe = (2 * self.eps_parallel + self.eps_perp) * (self.av)
# Kale's Thesis, Eq.(2.16)
self.Qe = (- self.b * (self.c11 + 2 * self.c12) / self.c11 *
self.eps_parallel)
# corrections to the method used to calculate band edges,
# thanks to Yu Song
# conduction band edge at different point, Eq.(2.7)
# Varsh correction is added later
self.EcG = self.VBO + self.EgG + self.Pec - bandBaseln
# band edge at L and X?
# only used in diagram..
self.EcL = (self.VBO + self.EgL +
(2 * self.eps_parallel + self.eps_perp) *
(self.acL + self.av) - bandBaseln)
self.EcX = (self.VBO + self.EgX +
(2 * self.eps_parallel + self.eps_perp) *
(self.acX + self.av) +
2 / 3 * self.XiX * (self.eps_perp - self.eps_parallel) -
bandBaseln)
self.ESO = sqrt(9 * self.Qe**2 + 2 * self.Qe * self.DSO + self.DSO**2)
self.EgLH = (self.EgG + self.Pec + self.Pe -
1 / 2 * (self.Qe - self.DSO + self.ESO))
self.EgSO = (self.EgG + self.Pec + self.Pe -
1 / 2 * (self.Qe - self.DSO - self.ESO))
# Varsh correction comes here
# temperature correction to conduction band edge, Eq.(2.10) in Kale's
self.Varsh = (- self.alG * cst.Temperature**2 /
(cst.Temperature + self.beG))
# the Varsh correction should be part conduction band, part valence
# 1st MAJOR assumption:
# Varshney contribution to band edge is in proportion to percent
# of band offset
# 2nd major assumption:
# Temperature affects sattelite valleys in the same way it does
# the Gamma valley
# ??But why is it working on offset between wells and barriers, not
# within the same material? This is not reasonable! (TODO)
# [?] Edward, T. Yu, James O. McCaldin, and Thomas C. McGill. "Band
# offsets in semiconductor heterojunctions." Solid state physics.
# Vol. 46. Academic Press, 1992. 1-146.
barrs = np.array([1, 3, 5, 7])
wells = np.array([0, 2, 4, 6])
CBOffset = self.EcG[barrs] - self.EcG[wells]
VBOffset = (self.EcG[barrs] - self.EgLH[barrs]) \
- (self.EcG[wells] - self.EgLH[wells])
percentCB = CBOffset / (CBOffset + VBOffset)
percentCB = np.column_stack([percentCB, percentCB]).flatten()
# applies percent CV to both well and barrier slots
self.EvLH = self.EcG - self.EgLH - ((1 - percentCB) * self.Varsh)
self.EvSO = self.EcG - self.EgSO - ((1 - percentCB) * self.Varsh)
self.EcG += percentCB * self.Varsh
self.EcL += percentCB * self.Varsh
self.EcX += percentCB * self.Varsh
# Eq.(2.20) in Kale's, with Eq=0. Note that E(C-SO) = EgSO = ESO+EgLH
self.me = 1 / ((1 + 2 * self.F) + self.Ep / self.EgLH * (
self.EgLH + 2 / 3 * self.ESO) / (self.EgLH + self.ESO))
def eff_mass(self, E):
"""Calculate effective mass according to energy E,
according to Eq.(2.20) in Kale's thesis
"""
# xMcE = self.xMc * (1 - (self.xVc - E) / self.xEg)
xMcE = 1 / (1 + 2 * self.xF + self.xEp / 3 * (
2 / ((E - self.xVc) + self.xEg) + 1 / (
(E - self.xVc) + self.xEg + self.xESO)))
return xMcE
def solve_psi(self):
""" solve eigen modes
OUTPUT: (doesn't return, but update member variables
self.EigenE is the eignenergy of the layer structure
self.xPointsPost[x] is a shorter version of self.xPoints
self.xyPsi[x, n] is the wave function at position
self.xPointsPost[x] corresiponding to the
eigenenergy EigenE[n], and without solutions near zero
self.xyPsiPsi[x, n] is the scaled norm of xyPsi
-- and the above two also cut long zero heads and tials --
-- for better plot --
self.xyPsiPsi2[x, n] is a more precise version corresponding to
position self.xPoints[x]
TODO: try matrix eigen solver?
"""
Epoints = np.arange(min(self.xVc),
max(self.xVc - 115 * self.EField * 1e-5),
self.vertRes / 1000)
xMcE = np.zeros(self.xPoints.shape)
xPsi = np.zeros(self.xPoints.shape)
psiEnd = np.zeros(Epoints.size)
# TODO: add adaptive spacing for Eq
if __USE_CLIB__:
# Call C function to get boundary dependence of energy EPoints[n],
# the return value is psiEnd[n]
# for n with psiEnd[n]=0, EPoints[n] is eigenenergy
cQ.psiFnEnd(Epoints.ctypes.data_as(ct.c_void_p),
int(Epoints.size), int(xPsi.size),
ct.c_double(self.xres), ct.c_double(self.EField),
self.xVc.ctypes.data_as(ct.c_void_p),
self.xEg.ctypes.data_as(ct.c_void_p),
self.xF.ctypes.data_as(ct.c_void_p),
self.xEp.ctypes.data_as(ct.c_void_p),
self.xESO.ctypes.data_as(ct.c_void_p),
self.xMc.ctypes.data_as(ct.c_void_p),
xMcE.ctypes.data_as(ct.c_void_p),
xPsi.ctypes.data_as(ct.c_void_p),
psiEnd.ctypes.data_as(ct.c_void_p))
else:
for p, Eq in enumerate(Epoints):
xMcE = m0 * self.eff_mass(Eq)
xMcE[0:-1] = 0.5 * (xMcE[0:-1] + xMcE[1:])
xPsi[0] = 0
xPsi[1] = 1
for q in xrange(1, xPsi.size - 1):
xPsi[q + 1] = xMcE[q] * (
(2 * (self.xres * ANG / hbar)**2 *
(self.xVc[q] - Eq) * e0 + 1 / xMcE[q] +
1 / xMcE[q - 1]) * xPsi[q] -
xPsi[q - 1] / xMcE[q - 1])
psiEnd[p] = xPsi[-1]
if __LOG__:
global logcount
with file("EpointsLog%d.pkl" % logcount, 'w') as logfile:
pickle.dump((Epoints, psiEnd), logfile)
logcount += 1
print "log saved for Epoints and psiEnd (%d)" % logcount
# TODO: maybe improved
tck = interpolate.splrep(Epoints, psiEnd)
self.EigenE = interpolate.sproot(tck, mest=len(Epoints))
if __MORE_INTERPOLATION__:
# Near the above approximation result,
# try to get a more precise result
xnear = np.empty(3 * len(self.EigenE))
fxnear = np.empty(3 * len(self.EigenE))
for q in xrange(self.EigenE.size):
# 100000 is an estimate for the precision of above
# approximation
# TODO: change the three calls of psiFn to loop
approxwidth = self.vertRes / 100000
xnear[3 * q: 3 * q + 3] = self.EigenE[q] + approxwidth * \
np.arange(-1, 2)
for n, Eq in enumerate(xnear):
if __USE_CLIB__:
cQ.psiFn(ct.c_double(Eq), int(1), int(xPsi.size),
ct.c_double(self.xres),
self.xVc.ctypes.data_as(ct.c_void_p),
self.xEg.ctypes.data_as(ct.c_void_p),
self.xF.ctypes.data_as(ct.c_void_p),
self.xEp.ctypes.data_as(ct.c_void_p),
self.xESO.ctypes.data_as(ct.c_void_p),
self.xMc.ctypes.data_as(ct.c_void_p),
xMcE.ctypes.data_as(ct.c_void_p),
xPsi.ctypes.data_as(ct.c_void_p))
else:
xMcE = self.eff_mass(Eq)
xMcE[1:] = (xMcE[:-1] + xMcE[1:]) / 2
xPsi[0] = 0
for x in range(1, xPsi.size):
# not tested.. C files are better
xPsi[x] = (xPsi[x] * (2 * (self.xres * ANG / hbar)**2 *
(self.xVc[x] - Eq) * e0 +
1 / xMcE[x] + 1 / xMcE[n - 1]) -
xPsi[x - 1] / xMcE[x - 1]) * xMcE[x]
fxnear[n] = xPsi[-1]
idxs = 3 * np.arange(len(self.EigenE)) + 1
if __USE_CLIB__:
cQ.inv_quadratic_interp(
xnear.ctypes.data_as(ct.c_void_p),
fxnear.ctypes.data_as(ct.c_void_p),
idxs.ctypes.data_as(ct.POINTER(ct.c_int)),
int(idxs.size), self.EigenE.ctypes.data_as(ct.c_void_p))
else:
for q, idx in enumerate(idxs):
# do quadratic interpolation
x0 = xnear[idx - 1]
fx0 = fxnear[idx - 1]
x1 = xnear[idx]
fx1 = fxnear[idx]
x2 = xnear[idx + 1]
fx2 = fxnear[idx + 1]
# inverse quadratic interpolation
x3 = (x0 * fx1 * fx2 / (fx0 - fx1) / (fx0 - fx2) +
x1 * fx0 * fx2 / (fx1 - fx0) / (fx1 - fx2) +
x2 * fx0 * fx1 / (fx2 - fx0) / (fx2 - fx1))
self.EigenE[q] = x3
# make array for Psi and fill it in
if __USE_CLIB__:
# with eigenenregy EigenE, here call C function to get wave
# function
self.xyPsi = np.zeros(self.xPoints.size * self.EigenE.size)
cQ.psiFill(int(xPsi.size), ct.c_double(self.xres),
int(self.EigenE.size),
self.EigenE.ctypes.data_as(ct.c_void_p),
self.xVc.ctypes.data_as(ct.c_void_p),
self.xEg.ctypes.data_as(ct.c_void_p),
self.xF.ctypes.data_as(ct.c_void_p),
self.xEp.ctypes.data_as(ct.c_void_p),
self.xESO.ctypes.data_as(ct.c_void_p),
self.xMc.ctypes.data_as(ct.c_void_p),
xMcE.ctypes.data_as(ct.c_void_p),
self.xyPsi.ctypes.data_as(ct.c_void_p))
# self.xyPsi.resize(a.xPoints.size, a.EigenE.size)
self.xyPsi = self.xyPsi.reshape(self.xPoints.size,
self.EigenE.size, order='F')
else:
self.xyPsi = np.zeros((self.xPoints.size, self.EigenE.size))
for p, Eq in enumerate(self.EigenE):
xMcE = m0 * self.eff_mass(self.EigenE)
xMcE[0:-1] = 0.5 * (xMcE[0:-1] + xMcE[1:])
xPsi[1] = 1
for q in xrange(2, xPsi.size - 1):
xPsi[q + 1] = ((2 * (self.xres * ANG / hbar)**2 *
(self.xVc[q] - Eq) * e0 + 1 / xMcE[q] +
1 / xMcE[q - 1]) * xPsi[q] -
xPsi[q - 1] / xMcE[q - 1]) * xMcE[q]
psiInt = np.sum(xPsi**2 * (
1 + (Eq - self.xVc) / (Eq - self.xVc + self.xEg)))
A = 1 / sqrt(self.xres * ANG * psiInt)
self.xyPsi[:, p] = A * xPsi
# remove states that come from oscillating end points
# TODO: change to remove non-bounded states, with user options
# wf with non-negeligiable amplitudes higher than barrier
# should be removed
if True:
# looks like we should change -1 to -2 (following)???
# psiEnd = self.xyPsi[-1,:]
# idxs = abs(psiEnd)<10
# idxs = np.nonzero(abs(psiEnd)<10)[0]
psiEnd = self.xyPsi[-2, :]
idxs = np.abs(psiEnd) < 200 / self.xres
# 200 depends on how precise we want about eigenenergy solver
# (TODO: more analysis and test about this value
self.EigenE = self.EigenE[idxs]
self.xyPsi = self.xyPsi[:, idxs]
# remove states that come from oscillating end points
# This doesn't really help for free levels, but accidentally removes
# the duplicate states on 0th layer.. (TODO)
psiEnd = self.xyPsi[-1, :]
idxs = np.nonzero(abs(psiEnd) < 10)[0]
self.EigenE = self.EigenE[idxs]
self.xyPsi = self.xyPsi[:, idxs]
# 4.5e-10 scales size of wavefunctions, arbitrary for nice plots
self.xyPsiPsi = self.xyPsi * self.xyPsi * wf_scale
# remove states that are smaller than minimum height (remove zero
# solutions?)-test case not showing any effect
# addresses states high above band edge
# 0.014 is arbitrary; changes if 4.5e-10 changes
# idxs = np.nonzero(self.xyPsiPsi.max(0) >
# wf_scale * wf_min_height)[0]
idxs = self.xyPsiPsi.max(0) > wf_scale * wf_min_height
self.EigenE = self.EigenE[idxs]
self.xyPsi = self.xyPsi[:, idxs]
self.xyPsiPsi = self.xyPsiPsi[:, idxs]
# self.xyPsiPsi2 = copy.deepcopy(self.xyPsiPsi)
self.xyPsiPlot = self.xyPsi * psi_scale
# implement pretty plot:
# remove long zero head and tail of the wave functions
# test case shows on "solve whole"
for q in xrange(self.EigenE.size):
# 0.0005 is arbitrary
prettyIdxs = np.nonzero(
self.xyPsiPsi[:, q] >
wf_scale * pretty_plot_factor)[0]
self.xyPsiPsi[0:prettyIdxs[0], q] = np.NaN
self.xyPsiPlot[0:prettyIdxs[0], q] = np.NaN
self.xyPsiPsi[prettyIdxs[-1]:, q] = np.NaN
self.xyPsiPlot[prettyIdxs[-1]:, q] = np.NaN
# decimate plot points: seems for better time and memory performance?
idxs = np.arange(0, self.xPoints.size, int(
plot_decimate_factor / self.xres), dtype=int)
# self.xyPsiPsiDec = np.zeros((idxs.size, self.EigenE.size))
# for q in xrange(self.EigenE.size):
# self.xyPsiPsiDec[:, q] = self.xyPsiPsi[idxs, q]
# self.xyPsiPsi = self.xyPsiPsiDec
self.xyPsiPsi = self.xyPsiPsi[idxs, :]
self.xyPsiPlot = self.xyPsiPlot[idxs, :]
self.xPointsPost = self.xPoints[idxs]
def solve_psi_mtr(self):
# TODO: a matrix eigen solver
pass
def basisSolve(self):
""" solve basis for the QC device, with each basis being eigen mode of
a seperate part of the layer structure
OUTPUT:
dCL: a list, each element is a QCLayers class, with layer structure
limited within a seperate sigle active/injection area, and
layer structure in dCL also includes pedding at head/tail
with same material as the first/last layer and barrier type
"""
# self.basisInjectorAR is 0-to-1
# self.basisARInjector is 1-to-0
# find all 0-to-1 & 1-to-0 transition points
# (1 for active region, 0 for injection retion, and active regions are
# required to start and end by a well layer)
# where for all n,
# self.layerARs[zeroTOone[n]] = self.layerARs[oneTOzero[n]] = 0
# but self.layerARs[zeroTOone[n]+1] = self.layerARs[oneTOzero[n]-1] = 1
# TODO: try always at left
zeroTOone = []
oneTOzero = []
layerAR = np.insert(self.layerARs, 0, self.layerARs[-1])
for q in xrange(0, layerAR.size - 1):
if layerAR[q] == 0 and layerAR[q + 1] == 1:
zeroTOone.append(q - 1)
if layerAR[q] == 1 and layerAR[q + 1] == 0:
oneTOzero.append(q)
# for q in xrange(0, self.layerARs.size - 1):
# if self.layerARs[q] == 0 and self.layerARs[q + 1] == 1:
# zeroTOone.append(q)
# if self.layerARs[q] == 1 and self.layerARs[q + 1] == 0:
# oneTOzero.append(q + 1)
dividers = [0, self.layerARs.size - 1]
if self.basisInjectorAR:
dividers += zeroTOone
if self.basisARInjector:
dividers += oneTOzero
# This converts the list into a set, thereby removing duplicates,
# and then back into a list.
dividers = list(set(dividers))
dividers.sort()
print dividers
# this is dataClassesList.
# it holds all of the Data classes for each individual solve section
dCL = []
# for first period only
# this handles all of the solving
for n in range(len(dividers) - 1):
dCL.append(copy.deepcopy(self))
dCL[n].repeats = 1
# substitute proper layer characteristics into dCL[n], hear/tail
# padding
layer = range(dividers[n], dividers[n + 1] + 1)
dCL[n].layerWidth = np.concatenate(
([int(PAD_HEAD / self.xres)], self.layerWidth[layer],
[int(PAD_TAIL / self.xres)]))
dCL[n].layerBarriers = np.concatenate(
([1], self.layerBarriers[layer], [1]))
dCL[n].layerARs = np.concatenate(
([0], self.layerARs[layer], [0]))
dCL[n].layerMaterials = np.concatenate(
([self.layerMaterials[layer][0]], self.layerMaterials[layer],
[self.layerMaterials[layer][-1]]))
dCL[n].layerDopings = np.concatenate(
([0], self.layerDopings[layer], [0]))
dCL[n].layerDividers = np.concatenate(
([0], self.layerDividers[layer], [0]))
# update and solve
dCL[n].update_alloys()
dCL[n].update_strain()
dCL[n].populate_x()
dCL[n].populate_x_band()
dCL[n].solve_psi()
print dCL[n].xVc
# caculate offsets
dCL[n].widthOffset = self.xres * np.sum(
self.layerWidth[range(0, dividers[n])])
dCL[n].fieldOffset = (-(dCL[n].widthOffset - PAD_HEAD) * ANG *
dCL[n].EField * KVpCM)
# create dCL's and offsets for repeat periods
period = len(dCL)
counter = period
if self.repeats > 1:
for q in xrange(1, self.repeats):
for p in xrange(0, period):
dCL.append(copy.deepcopy(dCL[p]))
dCL[counter].widthOffset = self.xres * np.sum(
self.layerWidth[1:]) * q + dCL[p].widthOffset
dCL[counter].fieldOffset = (-(dCL[counter].widthOffset -
100) * ANG *
dCL[counter].EField * KVpCM)
counter += 1
return dCL
def convert_dCL_to_data(self, dCL):
""" post processng of dCL list
INPUT:
self: orginal QCLayers class
dCL: result of basisSolve(self)
OUPUT:
get wave functions (dCL[n].xyPsi) and eigenenrgies (dCL[n].EigenE)
in dCL and update them in self; format them in length compatibale
for self and update self.xyPsiPsi
self.moduleID: moduleID[n] is the label of the position area for
mode self.eigenE[n] and self.xyPsi[n]
"""
# count number of wavefunctions
numWFs = sum([dC.EigenE.size for dC in dCL])
self.xPointsPost = np.arange(-PAD_HEAD, self.xPoints[-1] + PAD_TAIL +
self.xres, self.xres)
self.xyPsi = np.zeros((self.xPointsPost.size, numWFs))
self.xyPsiPsi = np.NaN * np.zeros(self.xyPsi.shape)
self.EigenE = np.zeros(numWFs)
self.moduleID = np.zeros(numWFs, dtype=np.int8)
counter = 0
for n, dC in enumerate(dCL):
for q in xrange(dC.EigenE.size):
self.EigenE[counter] = dC.EigenE[q] + dC.fieldOffset
self.moduleID[counter] = n
wf = np.zeros(self.xPointsPost.size)
begin = int(dC.widthOffset / self.xres)
end = begin + dC.xyPsi[:, q].size
wf[begin:end] = dC.xyPsi[:, q]
self.xyPsi[:, counter] = wf
self.xyPsiPsi[:, counter] = wf**2 * wf_scale
counter += 1
# cut head and tial to promise the figure is in the right place?
head = int(PAD_HEAD / self.xres)
tail = -int(PAD_TAIL / self.xres)
self.xPointsPost = self.xPointsPost[head:tail]
self.xyPsi = self.xyPsi[head:tail]
self.xyPsiPsi = self.xyPsiPsi[head:tail]
self.xyPsiPlot = self.xyPsi * psi_scale
# implement pretty plot
# remove long zero head and tail of the wave functions
# TODO: improve to cut according to range of well
for q in xrange(self.EigenE.size):
prettyIdxs = np.nonzero(
self.xyPsiPsi[:, q] >
wf_scale * pretty_plot_factor)[0]
if prettyIdxs.size != 0:
self.xyPsiPsi[0:prettyIdxs[0], q] = np.NaN
self.xyPsiPlot[0:prettyIdxs[0], q] = np.NaN
self.xyPsiPsi[prettyIdxs[-1]:, q] = np.NaN
self.xyPsiPlot[prettyIdxs[-1]:, q] = np.NaN
else:
self.xyPsiPsi[:, q] = np.NaN
self.xyPsiPlot[:, q] = np.NaN
# sort by ascending energy
sortID = np.argsort(self.EigenE)
self.EigenE = self.EigenE[sortID]
self.xyPsi = self.xyPsi[:, sortID]
self.xyPsiPsi = self.xyPsiPsi[:, sortID]
self.xyPsiPlot = self.xyPsiPlot[:, sortID]
self.moduleID = self.moduleID[sortID]
# #decimate plot points
# idxs = np.arange(0,self.xPoints.size, int(
# plot_decimate_factor/self.xres), dtype=int)
# self.xyPsiPsiDec = np.zeros([idxs.size, self.EigenE.size])
# for q in xrange(self.EigenE.size):
# self.xyPsiPsiDec[:,q] = self.xyPsiPsi[idxs,q]
# self.xyPsiPsi = self.xyPsiPsiDec
# self.xPointsPost = self.xPoints[idxs]
def lo_transition_rate(self, upper, lower):
""" LO phonon scattering induced decay life time calculator
INPUT:
upper: the higher energy state index
lower: the lower energy state index
OUTPUT
T1 decay life time between upper and lower states induced by LO
phonon scattering
"""
if upper < lower:
upper, lower = lower, upper
psi_i = self.xyPsi[:, upper]
psi_j = self.xyPsi[:, lower]
E_i = self.EigenE[upper]
E_j = self.EigenE[lower]
# print "---debug---"
# print "E_i = %f; E_j=%f"%(E_i, E_j)
if E_i - E_j - self.hwLO[0] < 0:
# energy difference is smaller than a LO phonon
# LO phonon scatering doesn't happen
return INV_INF
# zero head and tail cut off
idxs_i = np.nonzero(
psi_i > wf_scale * phonon_integral_factor)[0]