Effects of Doping on Finite Wells

Attached above is a PDF of my final project for Quantum Physics III with Professor Maxim Sukharev. This is where my knowledge of Fortran90 comes from and is how I gained knowledge of Akima Splines, DVR functions, Split operator functions, and many more implementations of quantum physics within Fortran90 while using Atom.io as an interface and Matlab for the graphs.

Below is a copy paste of the code used. It is written in Fortran 90

implicit none
double precision, parameter :: pi=3.141592653589793
double complex, parameter :: Im=(0.0d0,1.0d0) !imaginary unit
double precision, parameter :: hbar=1.0 !atomic units
integer, parameter :: N=2049 !total number of points
double precision, parameter :: a=-300.0,b=300.0 !coordinate interval (must be greater than potental width!)
double precision, parameter :: mass=1.0 !1 in atomic units means it is electron
double precision, parameter :: dx=(b-a)/(N-1)
double precision x
double precision H(N-1,N-1) !DVR representation of the Hamiltonian
double precision T,V,potential
double precision En(N-1) !eigenenergies
double precision psi(N-1,N-1) !eignefunctions
double precision const,tmp,tmp1,tmp2
double precision fv1(N-1),fv2(N-1) !working arrays for rs subroutine
integer i,j,ierr,matzconst=(hbar*pi)**2/(4.0*mass*(b-a)**2)
H=0.0
do i=1,N-1
x=a+dx*float(i)
do j=1,N-1
if(i==j)then
T=const*((2.0*(N**2)+1.0)/3.0-1.0/sin(pi*float(i)/float(N))**2)
V=potential(x)
else
T=const*((-1)**(i-j))*(1.0/sin(pi*float(i-j)/(2.0*float(N)))**2- &
1.0/sin(pi*float(i+j)/(2.0*float(N)))**2)
V=0.0
endif
H(i,j)=T+V !Hamiltonian
enddo
enddo
!H is ready
!call subroutine RS to diagonalize H
matz=1 !set this to 0 if ONLY eigenenergies needed
call rs(N-1,N-1,H,En,matz,psi,fv1,fv2,ierr) !
write(*,*) ierr !if ierr=0 => calculations have no errors, if not something went wrong
open(file='energies_doped_-1.dat',unit=11)
do i=1,50 !print first 50 energies
write(11,*) i,En(i)
enddo
close(unit=11)
!normalize all eignefunctions properly
!note that for psi(i,j) i is coordinate and j is quantum number
do i=1,N-1
tmp=0.0
do j=1,N-1
tmp=tmp+abs(psi(j,i))**2
enddo
tmp=tmp*dx
psi(:,i)=psi(:,i)/sqrt(tmp) !now it is propely normalized
enddo
open(file='states_doped_-1.dat',unit=11)
do i=1,N-1
x=a+dx*float(i)
write(11,*) x,psi(i,1),psi(i,2) !first two states
enddo
close(unit=11)
open(file='states34_doped_-1.dat',unit=11)
do i=1,N-1
x=a+dx*float(i)
write(11,*) x,psi(i,3),psi(i,4) !second two states
enddo
close(unit=11)
end
double precision function potential(x)
implicit none
double precision, parameter :: distance=10.0d0 !distance between wells
double precision, parameter :: width=10.0d0
double precision, parameter :: V0=-1.0
double precision, parameter :: V1=-1.1
double precision x
! if((x<distance/2.0).and.(x>-distance/2.0))then !well 0
! potential=V0
if((x>distance/2.0).and.(x<(width+distance/2.0)))then !well 1
potential=V0
elseif((x>(width+3.0*distance/2.0)).and.(x<(2.0*width+3.0*distance/2.0)))then !well 2
potential=V0
elseif((x>(2.0*width+5.0*distance/2.0)).and.(x<(3.0*width+5.0*distance/2.0)))then !well 3
potential=V0
elseif((x>(3.0*width+7.0*distance/2.0)).and.(x<(4.0*width+7.0*distance/2.0)))then !well 4
potential=V0
elseif((x>(4.0*width+9.0*distance/2.0)).and.(x<(5.0*width+9.0*distance/2.0)))then !well 5
potential=V0
elseif((x>(5.0*width+11.0*distance/2.0)).and.(x<(6.0*width+11.0*distance/2.0)))then !well 6
potential=V0
elseif((x>(6.0*width+13.0*distance/2.0)).and.(x<(7.0*width+13.0*distance/2.0)))then !well 7
potential=V0
elseif((x>(7.0*width+15.0*distance/2.0)).and.(x<(8.0*width+15.0*distance/2.0)))then !well 8
potential=V0
! elseif((x>(8.0*width+17.0*distance/2.0)).and.(x<(9.0*width+17.0*distance/2.0)))then !well 9
! potential=V0
! elseif((x>(9.0*width+19.0*distance/2.0)).and.(x<(10.0*width+19.0*distance/2.0)))then !well 10
! potential=V0
elseif((x<-distance/2.0).and.(x>-(width+distance/2.0)))then !well -1
potential=V0
elseif((x<-(width+3.0*distance/2.0)).and.(x>-(2.0*width+3.0*distance/2.0)))then !well -2
potential=V0
elseif((x<-(2.0*width+5.0*distance/2.0)).and.(x>-(3.0*width+5.0*distance/2.0)))then !well -3
potential=V0
elseif((x<-(3.0*width+7.0*distance/2.0)).and.(x>-(4.0*width+7.0*distance/2.0)))then !well -4
potential=V0
elseif((x<-(4.0*width+9.0*distance/2.0)).and.(x>-(5.0*width+9.0*distance/2.0)))then !well -5
potential=V0
elseif((x<-(5.0*width+11.0*distance/2.0)).and.(x>-(6.0*width+11.0*distance/2.0)))then !well -6
potential=V0
elseif((x<-(6.0*width+13.0*distance/2.0)).and.(x>-(7.0*width+13.0*distance/2.0)))then !well -7
potential=V0
elseif((x<-(7.0*width+15.0*distance/2.0)).and.(x>-(8.0*width+15.0*distance/2.0)))then !well -8
potential=V0
! elseif((x<-(8.0*width+17.0*distance/2.0)).and.(x>-(9.0*width+17.0*distance/2.0)))then !well -9
! potential=V0
! elseif((x<-(9.0*width+19.0*distance/2.0)).and.(x>-(10.0*width+19.0*distance/2.0)))then !well -10
! potential=V0
else
potential=0.0
endif
end function potential!double precision function potential(x)
! implicit none
! double precision, parameter :: distance=10.0d0 !distance between wells
! double precision, parameter :: width1=10.0d0 !width of well 1
! double precision, parameter :: width2=10.0d0 !width of well 2
! double precision a_V1,b_V1 !first well
! double precision a_V2,b_V2 !second well
! double precision, parameter :: V01=-100 !depth of the first well (all in atomic units)
! double precision, parameter :: V02=-100 !depth of the second well (all in atomic units)
! double precision, parameter :: Vmiddle=0.0 !height of the barrier between wells
! double precision x
! a_V1=-distance/2.0-width1
! b_V1=-distance/2.0
! a_V2=distance/2.0
! b_V2=distance/2.0+width2
! if((x>=a_v1).and.(x<=b_v1))then
! potential=V01
! elseif((x>=a_v2).and.(x<=b_v2))then
! potential=V02
! elseif((x>b_v1).and.(x<b_v2))then
! potential=Vmiddle
! else
! potential=0.0
! endif
!end function potential
subroutine rs(nm,n,a,w,matz,z,fv1,fv2,ierr)
integer n,nm,ierr,matz
double precision a(nm,n),w(n),z(nm,n),fv1(n),fv2(n)
!
! this subroutine calls the recommended sequence of
! subroutines from the eigensystem subroutine package (eispack)
! to find the eigenvalues and eigenvectors (if desired)
! of a real symmetric matrix.
!
! on input
!
! nm must be set to the row dimension of the two-dimensional
! array parameters as declared in the calling program
! dimension statement.
!
! n is the order of the matrix a.
!
! a contains the real symmetric matrix.
!
! matz is an integer variable set equal to zero if
! only eigenvalues are desired. otherwise it is set to
! any non-zero integer for both eigenvalues and eigenvectors.
!
! on output
!
! w contains the eigenvalues in ascending order.
!
! z contains the eigenvectors if matz is not zero.
!
! ierr is an integer output variable set equal to an error
! completion code described in the documentation for tqlrat
! and tql2. the normal completion code is zero.
!
! fv1 and fv2 are temporary storage arrays.
!
! questions and comments should be directed to burton s. garbow,
! mathematics and computer science div, argonne national laboratory
!
! this version dated august 1983.
!
! ------------------------------------------------------------------
!
if (n .le. nm) go to 10
ierr = 10 * n
go to 50
!
10 if (matz .ne. 0) go to 20
! .......... find eigenvalues only ..........
call tred1(nm,n,a,w,fv1,fv2)
! tqlrat encounters catastrophic underflow on the Vax
! call tqlrat(n,w,fv2,ierr)
call tql1(n,w,fv1,ierr)
go to 50
! .......... find both eigenvalues and eigenvectors ..........
20 call tred2(nm,n,a,w,fv1,z)
call tql2(nm,n,w,fv1,z,ierr)
50 return
end
double precision function pythag(a,b)
double precision a,b
!
! finds dsqrt(a**2+b**2) without overflow or destructive underflow
!
double precision p,r,s,t,u
p = dmax1(dabs(a),dabs(b))
if (p .eq. 0.0d0) go to 20
r = (dmin1(dabs(a),dabs(b))/p)**2
10 continue
t = 4.0d0 + r
if (t .eq. 4.0d0) go to 20
s = r/t
u = 1.0d0 + 2.0d0*s
p = u*p
r = (s/u)**2 * r
go to 10
20 pythag = p
return
end
subroutine tql1(n,d,e,ierr)
integer i,j,l,m,n,ii,l1,l2,mml,ierr
double precision d(n),e(n)
double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag
ierr = 0
if (n .eq. 1) go to 1001
!
do 100 i = 2, n
100 e(i-1) = e(i)
!
f = 0.0d0
tst1 = 0.0d0
e(n) = 0.0d0
!
do 290 l = 1, n
j = 0
h = dabs(d(l)) + dabs(e(l))
if (tst1 .lt. h) tst1 = h
! .......... look for small sub-diagonal element ..........
do 110 m = l, n
tst2 = tst1 + dabs(e(m))
if (tst2 .eq. tst1) go to 120
! .......... e(n) is always zero, so there is no exit
! through the bottom of the loop ..........
110 continue
!
120 if (m .eq. l) go to 210
130 if (j .eq. 30) go to 1000
j = j + 1
! .......... form shift ..........
l1 = l + 1
l2 = l1 + 1
g = d(l)
p = (d(l1) - g) / (2.0d0 * e(l))
r = pythag(p,1.0d0)
d(l) = e(l) / (p + dsign(r,p))
d(l1) = e(l) * (p + dsign(r,p))
dl1 = d(l1)
h = g - d(l)
if (l2 .gt. n) go to 145
!
do 140 i = l2, n
140 d(i) = d(i) - h
!
145 f = f + h
! .......... ql transformation ..........
p = d(m)
c = 1.0d0
c2 = c
el1 = e(l1)
s = 0.0d0
mml = m - l
! .......... for i=m-1 step -1 until l do -- ..........
do 200 ii = 1, mml
c3 = c2
c2 = c
s2 = s
i = m - ii
g = c * e(i)
h = c * p
r = pythag(p,e(i))
e(i+1) = s * r
s = e(i) / r
c = p / r
p = c * d(i) - s * g
d(i+1) = h + s * (c * g + s * d(i))
200 continue
!
p = -s * s2 * c3 * el1 * e(l) / dl1
e(l) = s * p
d(l) = c * p
tst2 = tst1 + dabs(e(l))
if (tst2 .gt. tst1) go to 130
210 p = d(l) + f
! .......... order eigenvalues ..........
if (l .eq. 1) go to 250
! .......... for i=l step -1 until 2 do -- ..........
do 230 ii = 2, l
i = l + 2 - ii
if (p .ge. d(i-1)) go to 270
d(i) = d(i-1)
230 continue
!
250 i = 1
270 d(i) = p
290 continue
!
go to 1001
! .......... set error -- no convergence to an
! eigenvalue after 30 iterations ..........
1000 ierr = l
1001 return
end
subroutine tql2(nm,n,d,e,z,ierr)
integer i,j,k,l,m,n,ii,l1,l2,nm,mml,ierr
double precision d(n),e(n),z(nm,n)
double precision c,c2,c3,dl1,el1,f,g,h,p,r,s,s2,tst1,tst2,pythag
ierr = 0
if (n .eq. 1) go to 1001
do 100 i = 2, n
100 e(i-1) = e(i)
f = 0.0d0
tst1 = 0.0d0
e(n) = 0.0d0
do 240 l = 1, n
j = 0
h = dabs(d(l)) + dabs(e(l))
if (tst1 .lt. h) tst1 = h
do 110 m = l, n
tst2 = tst1 + dabs(e(m))
if (tst2 .eq. tst1) go to 120
110 continue
120 if (m .eq. l) go to 220
130 if (j .eq. 30) go to 1000
j = j + 1
l1 = l + 1
l2 = l1 + 1
g = d(l)
p = (d(l1) - g) / (2.0d0 * e(l))
r = pythag(p,1.0d0)
d(l) = e(l) / (p + dsign(r,p))
d(l1) = e(l) * (p + dsign(r,p))
dl1 = d(l1)
h = g - d(l)
if (l2 .gt. n) go to 145
do 140 i = l2, n
140 d(i) = d(i) - h
145 f = f + h
p = d(m)
c = 1.0d0
c2 = c
el1 = e(l1)
s = 0.0d0
mml = m - l
do 200 ii = 1, mml
c3 = c2
c2 = c
s2 = s
i = m - ii
g = c * e(i)
h = c * p
r = pythag(p,e(i))
e(i+1) = s * r
s = e(i) / r
c = p / r
p = c * d(i) - s * g
d(i+1) = h + s * (c * g + s * d(i))
do 180 k = 1, n
h = z(k,i+1)
z(k,i+1) = s * z(k,i) + c * h
z(k,i) = c * z(k,i) - s * h
180 continue
200 continue
p = -s * s2 * c3 * el1 * e(l) / dl1
e(l) = s * p
d(l) = c * p
tst2 = tst1 + dabs(e(l))
if (tst2 .gt. tst1) go to 130
220 d(l) = d(l) + f
240 continue
do 300 ii = 2, n
i = ii - 1
k = i
p = d(i)
do 260 j = ii, n
if (d(j) .ge. p) go to 260
k = j
p = d(j)
260 continue
if (k .eq. i) go to 300
d(k) = d(i)
d(i) = p
do 280 j = 1, n
p = z(j,i)
z(j,i) = z(j,k)
z(j,k) = p
280 continue
300 continue
go to 1001
1000 ierr = l
1001 return
end
subroutine tred1(nm,n,a,d,e,e2)
integer i,j,k,l,n,ii,nm,jp1
double precision a(nm,n),d(n),e(n),e2(n)
double precision f,g,h,scale
do 100 i = 1, n
d(i) = a(n,i)
a(n,i) = a(i,i)
100 continue
do 300 ii = 1, n
i = n + 1 - ii
l = i - 1
h = 0.0d0
scale = 0.0d0
if (l .lt. 1) go to 130
do 120 k = 1, l
120 scale = scale + dabs(d(k))
if (scale .ne. 0.0d0) go to 140
do 125 j = 1, l
d(j) = a(l,j)
a(l,j) = a(i,j)
a(i,j) = 0.0d0
125 continue
130 e(i) = 0.0d0
e2(i) = 0.0d0
go to 300
140 do 150 k = 1, l
d(k) = d(k) / scale
h = h + d(k) * d(k)
150 continue
e2(i) = scale * scale * h
f = d(l)
g = -dsign(dsqrt(h),f)
e(i) = scale * g
h = h - f * g
d(l) = f - g
if (l .eq. 1) go to 285
do 170 j = 1, l
170 e(j) = 0.0d0
do 240 j = 1, l
f = d(j)
g = e(j) + a(j,j) * f
jp1 = j + 1
if (l .lt. jp1) go to 220
do 200 k = jp1, l
g = g + a(k,j) * d(k)
e(k) = e(k) + a(k,j) * f
200 continue
220 e(j) = g
240 continue
f = 0.0d0
do 245 j = 1, l
e(j) = e(j) / h
f = f + e(j) * d(j)
245 continue
h = f / (h + h)
do 250 j = 1, l
250 e(j) = e(j) - h * d(j)
do 280 j = 1, l
f = d(j)
g = e(j)
do 260 k = j, l
260 a(k,j) = a(k,j) - f * e(k) - g * d(k)
280 continue
285 do 290 j = 1, l
f = d(j)
d(j) = a(l,j)
a(l,j) = a(i,j)
a(i,j) = f * scale
290 continue
300 continue
return
end
subroutine tred2(nm,n,a,d,e,z)
integer i,j,k,l,n,ii,nm,jp1
double precision a(nm,n),d(n),e(n),z(nm,n)
double precision f,g,h,hh,scale
do 100 i = 1, n
do 80 j = i, n
80 z(j,i) = a(j,i)
d(i) = a(n,i)
100 continue
if (n .eq. 1) go to 510
do 300 ii = 2, n
i = n + 2 - ii
l = i - 1
h = 0.0d0
scale = 0.0d0
if (l .lt. 2) go to 130
do 120 k = 1, l
120 scale = scale + dabs(d(k))
if (scale .ne. 0.0d0) go to 140
130 e(i) = d(l)
do 135 j = 1, l
d(j) = z(l,j)
z(i,j) = 0.0d0
z(j,i) = 0.0d0
135 continue
go to 290
140 do 150 k = 1, l
d(k) = d(k) / scale
h = h + d(k) * d(k)
150 continue
f = d(l)
g = -dsign(dsqrt(h),f)
e(i) = scale * g
h = h - f * g
d(l) = f - g
do 170 j = 1, l
170 e(j) = 0.0d0
do 240 j = 1, l
f = d(j)
z(j,i) = f
g = e(j) + z(j,j) * f
jp1 = j + 1
if (l .lt. jp1) go to 220
do 200 k = jp1, l
g = g + z(k,j) * d(k)
e(k) = e(k) + z(k,j) * f
200 continue
220 e(j) = g
240 continue
f = 0.0d0
do 245 j = 1, l
e(j) = e(j) / h
f = f + e(j) * d(j)
245 continue
hh = f / (h + h)
do 250 j = 1, l
250 e(j) = e(j) - hh * d(j)
do 280 j = 1, l
f = d(j)
g = e(j)
do 260 k = j, l
260 z(k,j) = z(k,j) - f * e(k) - g * d(k)
d(j) = z(l,j)
z(i,j) = 0.0d0
280 continue
290 d(i) = h
300 continue
do 500 i = 2, n
l = i - 1
z(n,l) = z(l,l)
z(l,l) = 1.0d0
h = d(i)
if (h .eq. 0.0d0) go to 380
do 330 k = 1, l
330 d(k) = z(k,i) / h
do 360 j = 1, l
g = 0.0d0
do 340 k = 1, l
340 g = g + z(k,i) * z(k,j)
do 360 k = 1, l
z(k,j) = z(k,j) - g * d(k)
360 continue
380 do 400 k = 1, l
400 z(k,i) = 0.0d0
500 continue
510 do 520 i = 1, n
d(i) = z(n,i)
z(n,i) = 0.0d0
520 continue
z(n,n) = 1.0d0
e(1) = 0.0d0
return
end