aestimosolver / aestimo

Aestimo 1D Schrödinger-Poisson Solver
https://aestimosolver.github.io
GNU General Public License v3.0
50 stars 24 forks source link

InGaAs quantum well problem #32

Closed thomasrockett closed 2 years ago

thomasrockett commented 2 years ago

Hi I am a new user to both python and aestimo.

I am getting some strange results when trying to calculate the bound state energies in an InGaAs quantum well. As a sanity check I calculated the energy levels in an InGaAs well (with 0% indium content) and a GaAs well at the same time. These should give the same result, however the InGaAs bound state energies are hovering somewhere in the forbidden gap. Is this a bug, or due to an incorrect value in the database for InGaAs? Apologies if there is a really obvious solution as I haven't used python before.

I have a few other questions which I might as well include here:

Many thanks for viewing this post, this software looks extremely useful and I would be very happy to cite the recent paper in my own papers. Tom structure

sblisesivdin commented 2 years ago

You may use Quantum_Regions to calculate each QW as independent of each other. Otherwise, the whole system will be solved like a single coupled system.

One example: https://github.com/aestimosolver/aestimo/blob/11d7d067be647734cf4abced9c563b5adc32790f/examples/sample_qw_barrierdope_p_algan_gan_heterojunction.py

You can view other examples for different use cases.

thomasrockett commented 2 years ago

Hi Sefer

Thanks for your fast response to my problem and your patience with me. I am really sorry but I am still unable to get a realistic result from my simulation, even when I simulate the 2 wells separately. The bound states for the InGaAs (0% In) well are inside the band gap. I also noticed that the h2 wavefunction in the GaAs well looks symmetric, this should be an antisymmetric wavefunction as far as I know.

I have attached my code to this post, would it be possible for you to take a look and see if there is anything obvious I am doing wrong? I copied nearly all the code from the sample double qw example file.

Many thanks, Tom

Figure_2

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------
# Input File Description:  Double Quantum well, comparing InGaAs (x=0%) to GaAs quantum wells
# ------------------------------------------------------------------------
# ----------------
# GENERAL SETTINGS
# ----------------

# TEMPERATURE
T = 300.0 #Kelvin

# COMPUTATIONAL SCHEME
# 0: Schrodinger
# 1: Schrodinger + nonparabolicity
# 2: Schrodinger-Poisson
# 3: Schrodinger-Poisson with nonparabolicity
# 4: Schrodinger-Exchange interaction
# 5: Schrodinger-Poisson + Exchange interaction
# 6: Schrodinger-Poisson + Exchange interaction with nonparabolicity
# 7: Schrodinger-Poisson-Drift_Diffusion (Schrodinger solved with poisson then  poisson and DD)
# 8: Schrodinger-Poisson-Drift_Diffusion (Schrodinger solved with poisson and DD)
# 9: Schrodinger-Poisson-Drift_Diffusion (Schrodinger solved with poisson and DD) using Gummel & Newton map
computation_scheme = 2
# Non-parabolic effective mass function
# 0: no energy dependence
# 1: Nelson's effective 2-band model
# 2: k.p model from Vurgaftman's 2001 paper
meff_method = 2

# Non-parabolic Dispersion Calculations for Fermi-Dirac
fermi_np_scheme = True

# QUANTUM
# Total subband number to be calculated for electrons
subnumber_e = 2
subnumber_h = 2
# APPLIED ELECTRIC FIELD
Fapplied = 0.0 # (V/m)
vmax= 1.7
vmin= 0.0
Each_Step=0.05
mat_type='Zincblende'
# --------------------------------
# REGIONAL SETTINGS FOR SIMULATION
# --------------------------------

# GRID
# For 1D, z-axis is choosen
gridfactor = 0.5 #nm
maxgridpoints = 200000 #for controlling the size
# REGIONS
# Region input is a two-dimensional list input.
# An example:
# Si p-n diode. Firstly lets picturize the regional input.
#         | Thickness (nm)  | Material | Alloy fraction | Doping(cm^-3) | n or p type |
# Layer 0 |       250.0     |   Si     |      0         |     1e16      |     n       |
# Layer 1 |       250.0     |   Si     |      0         |     1e16      |     p       |
#
# To input this list in Gallium, we use lists as:

material =[[ 50, 'AlGaAs', 0.3, 0.0,0.0, 'i','NA'],
           [ 10, 'InGaAs', 0.0, 0.0,0.0, 'i','NA'],
           [ 50, 'AlGaAs',0.3, 0.0, 0.0, 'i','NA'],
           [ 10, 'GaAs', 0, 0.0,0.0, 'i','NA'],
           [50, 'AlGaAs', 0.3, 0.0, 0.0, 'i', 'NA']]

#----------------------------------------
#Doping profiles based on the LSS theory (ion implantation).
import numpy as np
x_max = sum([layer[0] for layer in material])
def round2int(x):
    return int(x+0.5)
n_max=round2int(x_max/gridfactor)
#----------------------------------------
dop_n=np.zeros(n_max)
dop_p=np.zeros(n_max)
dop_profile=np.zeros(n_max)
xaxis = np.arange(0,n_max)*gridfactor#[nm]
Q_n=2e12#implant dose [1/cm2]
Rp_n=86#projected range Rp [nm]
Delta_Rp_n=44#projected straggle Delta Rp [nm]
Q_p=1e11#implant dose [1/cm2]
Rp_p=75#projected range Rp [nm]
Delta_Rp_p=20#projected straggle Delta Rp [nm]
from math import sqrt, exp
def Lss_profile_dop(x,Q,Delta_Rp,Rp):
    return Q/(sqrt(2*np.pi)*Delta_Rp*1e-7)*exp(-(x-Rp)**2/(2*Delta_Rp**2))
def Lss_profile_dop_diff(x,Q,Delta_Rp,Rp):
    return Q/(2*sqrt(np.pi)*Delta_Rp*1e-7)*exp(-(x-Rp)**2/(4*Delta_Rp**2))
for i in range(n_max):
    dop_n[i]=Lss_profile_dop(xaxis[n_max-1-i],Q_n,Delta_Rp_n,Rp_n)*1e6
    dop_p[i]=-Lss_profile_dop(xaxis[n_max-1-i],Q_p,Delta_Rp_p,Rp_p)*1e6
    #dop_profile[i]=dop_n[i]+dop_p[i]
#----------------------------------------
#Quantum_Regions=False
#Quantum_Regions_boundary=np.zeros((2,2))

Quantum_Regions=True
#--------------------------------
Quantum_Regions_boundary=np.zeros((2,2))
Quantum_Regions_boundary[0,0]=30
Quantum_Regions_boundary[0,1]=80

Quantum_Regions_boundary[1,0]=90
Quantum_Regions_boundary[1,1]=140

#----------------------------------------
surface=np.zeros(2)
#surface[0]=-0.6
#----------------------------------------
if __name__ == "__main__": #this code allows you to run the input file directly
    input_obj = vars()
    import aestimo_eh
    aestimo_eh.run_aestimo(input_obj)

# Output Files
# ------------
output_directory = "output_"+inputfilename
parameters = True
electricfield_out = True
potential_out = True
sigma_out = True
probability_out = True
states_out = True
Drift_Diffusion_out=True
h-hebal commented 2 years ago

Dear Thomas,If you just write GaAs instead of InGaAs the simulation works fine, you have to check the parameters of InGaAs alloy in database, the second ratio  is for "y" in quaternary materials.RegardsHamza

Envoyé depuis Yahoo Mail pour Android

Le mar., janv. 11, 2022 à 11:17, @.***> a écrit:

Hi Sefer

Thanks for your fast response to my problem and your patience with me. I am really sorry but I am still unable to get a realistic result from my simulation, even when I simulate the 2 wells separately. The bound states for the InGaAs (0% In) well are inside the band gap. I also noticed that the h2 wavefunction in the GaAs well looks symmetric, this should be an antisymmetric wavefunction as far as I know.

I have attached my code to this post, would it be possible for you to take a look and see if there is anything obvious I am doing wrong? I copied nearly all the code from the sample double qw example file.

Many thanks, Tom

!/usr/bin/env python

-- coding: utf-8 --

------------------------------------------------------------------------

Input File Description: Double Quantum well, comparing InGaAs (x=0%) to GaAs quantum wells

------------------------------------------------------------------------

----------------

GENERAL SETTINGS

----------------

TEMPERATURE

T = 300.0 #Kelvin

COMPUTATIONAL SCHEME

0: Schrodinger

1: Schrodinger + nonparabolicity

2: Schrodinger-Poisson

3: Schrodinger-Poisson with nonparabolicity

4: Schrodinger-Exchange interaction

5: Schrodinger-Poisson + Exchange interaction

6: Schrodinger-Poisson + Exchange interaction with nonparabolicity

7: Schrodinger-Poisson-Drift_Diffusion (Schrodinger solved with poisson then poisson and DD)

8: Schrodinger-Poisson-Drift_Diffusion (Schrodinger solved with poisson and DD)

9: Schrodinger-Poisson-Drift_Diffusion (Schrodinger solved with poisson and DD) using Gummel & Newton map

computation_scheme = 2

Non-parabolic effective mass function

0: no energy dependence

1: Nelson's effective 2-band model

2: k.p model from Vurgaftman's 2001 paper

meff_method = 2

Non-parabolic Dispersion Calculations for Fermi-Dirac

fermi_np_scheme = True

QUANTUM

Total subband number to be calculated for electrons

subnumber_e = 2 subnumber_h = 2

APPLIED ELECTRIC FIELD

Fapplied = 0.0 # (V/m) vmax= 1.7 vmin= 0.0 Each_Step=0.05 mat_type='Zincblende'

--------------------------------

REGIONAL SETTINGS FOR SIMULATION

--------------------------------

GRID

For 1D, z-axis is choosen

gridfactor = 0.5 #nm maxgridpoints = 200000 #for controlling the size

REGIONS

Region input is a two-dimensional list input.

An example:

Si p-n diode. Firstly lets picturize the regional input.

| Thickness (nm) | Material | Alloy fraction | Doping(cm^-3) | n or p type |

Layer 0 | 250.0 | Si | 0 | 1e16 | n |

Layer 1 | 250.0 | Si | 0 | 1e16 | p |

#

To input this list in Gallium, we use lists as:

material =[[ 50, 'AlGaAs', 0.3, 0.0,0.0, 'i','NA'], [ 10, 'InGaAs', 0.0, 0.0,0.0, 'i','NA'], [ 50, 'AlGaAs',0.3, 0.0, 0.0, 'i','NA'], [ 10, 'GaAs', 0, 0.0,0.0, 'i','NA'], [50, 'AlGaAs', 0.3, 0.0, 0.0, 'i', 'NA']]

----------------------------------------

Doping profiles based on the LSS theory (ion implantation).

import numpy as np x_max = sum([layer[0] for layer in material]) def round2int(x): return int(x+0.5) n_max=round2int(x_max/gridfactor)

----------------------------------------

dop_n=np.zeros(n_max) dop_p=np.zeros(n_max) dop_profile=np.zeros(n_max) xaxis = np.arange(0,n_max)gridfactor#[nm] Q_n=2e12#implant dose [1/cm2] Rp_n=86#projected range Rp [nm] Delta_Rp_n=44#projected straggle Delta Rp [nm] Q_p=1e11#implant dose [1/cm2] Rp_p=75#projected range Rp [nm] Delta_Rp_p=20#projected straggle Delta Rp [nm] from math import sqrt, exp def Lss_profile_dop(x,Q,Delta_Rp,Rp): return Q/(sqrt(2np.pi)Delta_Rp1e-7)*exp(-(x-Rp)2/(2*Delta_Rp2)) def Lss_profile_dop_diff(x,Q,Delta_Rp,Rp): return Q/(2sqrt(np.pi)Delta_Rp1e-7)exp(-(x-Rp)2/(4*Delta_Rp2)) for i in range(n_max): dop_n[i]=Lss_profile_dop(xaxis[n_max-1-i],Q_n,Delta_Rp_n,Rp_n)1e6 dop_p[i]=-Lss_profile_dop(xaxis[n_max-1-i],Q_p,Delta_Rp_p,Rp_p)1e6

dop_profile[i]=dop_n[i]+dop_p[i]

----------------------------------------

Quantum_Regions=False

Quantum_Regions_boundary=np.zeros((2,2))

Quantum_Regions=True

--------------------------------

Quantum_Regions_boundary=np.zeros((2,2)) Quantum_Regions_boundary[0,0]=30 Quantum_Regions_boundary[0,1]=80

Quantum_Regions_boundary[1,0]=90 Quantum_Regions_boundary[1,1]=140

----------------------------------------

surface=np.zeros(2)

surface[0]=-0.6

----------------------------------------

if name == "main": #this code allows you to run the input file directly input_obj = vars() import aestimo_eh aestimo_eh.run_aestimo(input_obj)

Output Files

------------

outputdirectory = "output"+inputfilename parameters = True electricfield_out = True potential_out = True sigma_out = True probability_out = True states_out = True Drift_Diffusion_out=True

— Reply to this email directly, view it on GitHub, or unsubscribe. Triage notifications on the go with GitHub Mobile for iOS or Android. You are receiving this because you are subscribed to this thread.Message ID: @.***>

thomasrockett commented 2 years ago

Hi Hamza

Thank you for your reply. I will try to fix the database values for InGaAs but my background is in MBE growth so I am not sure what some of these parameters represent. Do you have any tips as to which parameters could be responsible for these strange results I am seeing? It seems that InGaAs and GaAsP both give strange results, but AlGaAs seems to be working, see picture.

Many thanks, Tom

Figure_3

thomasrockett commented 2 years ago

Update: I have found the problem with the database. I went through the InGaAs and InAs parameters one by one as you suggested, and found that the lattice constant a0 is responsible. Setting the InAs lattice constant to 5.66 Angstroms in the database fixes the problem, and I obtain the expected result! This seems slightly strange to me, so I will defer to your expertise on whether this constitutes a bug. GaAsP also seems to work if I change the GaP lattice constant to 5.66 Angstroms.

h-hebal commented 2 years ago

Dear Thomas, thank you for your fellowup, your are right after checking we found a misprint in the code regarding the linear estimation of the lattice constant in ternary zinkblende alloys and a missing variable : In aestimo_eh.py:

1-first correction if mat_crys_strc == "Zincblende": . . . a0[startindex:finishindex] = (x mat1["a0"] + (1 - x) mat2["a0"] before was: a0[startindex:finishindex] = ( (1 - x) mat1["a0"] +x mat2["a0"] 2- second correction: if mat_crys_strc == "Zincblende": a0_sub[startindex:finishindex] = matprops["a0_sub"] 1e-10 before was: a0_sub[startindex:finishindex] = matprops["a0"] 1e-10 Best regards Hamza

h-hebal commented 2 years ago

thomasrockett I think the issue is fixed therefore i am going to close it.