C ====================================================================== C PROGRAM 1DSP.f C 1D Eigenvalue solver that solves for the eigenvalues C and eigenvectors for arbitrary potential C Last modified: September 20th, 1998 C ====================================================================== IMPLICIT REAL*8(A-H, O-Z) PARAMETER(npts = 1000, n_sub_max = 20) REAL*8 q, h, hb, k_b, eps_o, mo DATA q, h, hb, k_b, eps_o, mo 1 / 1.602D-19, 6.625D-34, 1.0544D-34, 2 1.38D-23, 8.86D-12, 9.1D-31/ REAL*8 z_grid(0:npts) REAL*8 V_vec(0:npts) REAL*8 psi_vec(0:npts,n_sub_max) REAL*8 psi(0:npts) REAL*8 e_sub(n_sub_max) REAL*8 energy(0:npts,n_sub_max) C----- Description of important parameters: C npts => total number of points in the mesh C n_sub_max => maximum number of eigenvalues which will be searched C z_grid => 1D vector that contains the mesh C V_vec => Vector for the confining potential at the mesh points in eV C psi_vec => vector in which we store the wavefunctions C e_sub => vector with subband eneries C rmt => effective mass C n_sub => user determined number of subbands that will be calculated C i_fix => number of subbands found C subroutine RNUMEROV: numerical integration of the 1D Schrodinger equation C subroutine SCHRED: eigenvalue surch by using the bisection method C----- Set up the grid and the confining potential: ------------------------ print*,'Define well:1=square,2=parabolic,3=triangular,4=V-shaped' read*,flag n_sub = n_sub_max print*,'Effective mass' read*,rmt rmt = rmt*mo rmt = 6.D-32 dz = 2.D-10 zz = -float(npts)*dz/2 open(unit = 1, file = 'wavefunctions', status = 'unknown') open(unit = 2, file = 'energies', status = 'unknown') open(unit = 3, file = 'separate_wavefunctions', status = 'unknown') if(flag.eq.1)then print*,'well width [nm]' read*,well_length well_length = well_length*1D-9 print*,'well depth [eV]' read*,well_height endif if(flag.eq.2)then print*,'oscillator energy [eV]' read*,Energy_osc endif if(flag.eq.3)then print*,'Electric field [V/m]' read*,Electric_field endif if(flag.eq.4)then print*,'Electric field [V/m]' read*,Electric_field endif DO J = 0, npts z_grid(j) = dz zz = zz + dz if(flag.eq.1)then if(dabs(zz).gt.(well_length/2.D0))then V_vec(j) = well_height else V_vec(j) = 0. endif endif if(flag.eq.2)then term = 0.5*rmt*Energy_osc*Energy_osc/hb*zz/hb*zz*q V_vec(j) = term endif if(flag.eq.3)then V_vec(j) = Electric_field*float(j)*dz endif if(flag.eq.4)then V_vec(j) = Electric_field*float(j)*abs(zz) endif ENDDO C----- Calculate the subband structure: ------------------------------------- if(n_sub.ne.0)then CALL SCHRED(V_vec, rmt, z_grid, psi_vec, e_sub, 1 n_sub, i_fix, flag_end) print*,i_fix endif zz = -float(npts)*dz/2 do j = 0, npts write(1,2)zz,V_vec(j),((psi_vec(j,k)+e_sub(k)),k=1,i_fix) 2 format(22(E10.3,2X)) zz = zz + dz enddo do i = 1,i_fix do j = 0,npts term = psi_vec(j,i)**2 if(term.gt.(0.1))then energy(j,i)=e_sub(i) endif enddo enddo zz = -float(npts)*dz/2 do j = 0, npts write(2,3)zz,V_vec(j),(energy(j,k),k=1,i_fix) 3 format(22(E10.3,2X)) zz = zz + dz enddo zz = -float(npts)*dz/2 do j = 0, npts write(3,4)zz,V_vec(j),(psi_vec(j,k),k=1,i_fix) 4 format(22(E10.3,2X)) zz = zz + dz enddo END C ======================================================================= C SUBROUTINE FOR THE EVALUATION OF THE SUBBAND STRUCTURE BY USING THE C SHOOTING METHOD C The numerical integration of the Schroedinger equation is implemented C in the subroutine NUMEROV. C ======================================================================= SUBROUTINE SCHRED(V, md, z_grid, psi1, e_sub, 1 i_sub_max, i_fix, flag_end) IMPLICIT REAL*8(A-H, O-Z) PARAMETER(npts = 1000, n_sub_max = 20) REAL*8 q, h, hb, k_b, eps_o, mo DATA q, h, hb, k_b, eps_o, mo 1 / 1.602D-19, 6.625D-34, 1.0544D-34, 2 1.38D-23, 8.86D-12, 9.1D-31/ REAL*8 V(0:npts) REAL*8 z_grid(0:npts) REAL*8 e_sub(n_sub_max) REAL*8 psi(0:npts) REAL*8 psi1(0:npts,n_sub_max) REAL*8 psi1_max(n_sub_max), psi1_sum(n_sub_max) REAL*8 md REAL*8 fx(0:10000) logical I\$FLAG, FLAG_SCHRED dE_old = 0.001 dE_start = dE_old dE = dE_old tolf = 1.D-8 i_sub = 0 i_fix = 0 i_max = npts V_min = V(0) do i = 1, npts if(V(i).lt.V_min)V_min = V(i) enddo V_max = V(0) do i = 1, npts if(V(i).gt.V_max)V_max = V(i) enddo FLAG_SCHRED = .FALSE. e = V_min e = e + dE call rnumerov(V,md,e,z_grid,i_max,f,psi) fx(0) = f DO WHILE(.NOT.FLAG_SCHRED) I\$FLAG = .FALSE. k = 0 DO WHILE(.NOT.I\$FLAG) k = k + 1 if(k.ge.10000)then write(6,*)'There is a problem in the eigenvalue search!' write(6,*)'Try reducing the number of eigenvalues!' flag_end = 1. endif e = e + dE if(e.ge.V_max)then i_fix = i_sub FLAG_SCHRED = .TRUE. goto 21 endif call rnumerov(V,md,e,z_grid,i_max,f,psi) c print*,e,f fx(k) = f IF(fx(k-1)*fx(k).le.0)THEN e1 = e - dE I\$FLAG = .TRUE. ENDIF ENDDO rtbis = e1 I\$FLAG = .FALSE. DO WHILE(.NOT.I\$FLAG) dE = 0.5D0*dE e_mid = rtbis + dE call rnumerov(V,md,e_mid,z_grid,i_max,f,psi) fmid = f c print*,e_mid,fmid if(fx(0).lt.0)then if(fmid.lt.0.D0)rtbis = e_mid elseif(fx(0).gt.0)then if(fmid.gt.0.D0)rtbis = e_mid endif if(dE.lt.1.D-30.or.dabs(fmid).lt.tolf)I\$FLAG = .TRUE. ENDDO i_sub = i_sub + 1 e_sub(i_sub) = e_mid do i = 0, i_max if((psi(1).gt.0.).and.(mod(i_sub,2).eq.0))then psi1(i,i_sub) = - psi(i) else psi1(i,i_sub) = psi(i) endif enddo print*,'i_sub=',i_sub,' e=',e_mid,' error =',fmid if(i_sub.eq.i_sub_max)then i_fix = i_sub_max FLAG_SCHRED = .TRUE. endif if(i_sub.ge.2)then dE = (e_sub(i_sub)-e_sub(i_sub-1))/10.D0 else dE_old = 0.5*dE_old dE = dE_old endif e = e_mid + dE call rnumerov(V,md,e,z_grid,i_max,f,psi) fx(0) = f 21 continue ENDDO C-----Renormalize the evaluated wavefunctions to value 1 do j=1,i_fix psi1_max(j) = 0. do i = 0, i_max if(dabs(psi1(i,j)).ge.psi1_max(j)) 1 psi1_max(j)=dabs(psi1(i,j)) enddo if(psi1_max(j).ne.0.)then do i = 0, i_max psi1(i,j) = psi1(i,j)/psi1_max(j) enddo endif enddo C-----Calculate the square of the integral of the wavefunctions c do j = 1, i_fix c psi1_sum(j) = 0. c do i = 1, i_max c term = 0.5*(psi1(i-1,j) + psi1(i,j)) c psi1_sum(j)= psi1_sum(j) + term*term*z_grid(i-1) c enddo c enddo C-----Renormalize the wavefunctions to get an orthonormal set c do j = 1, i_fix c if(psi1_sum(j).ne.0.)then c do i = 0, i_max c psi1(i,j) = psi1(i,j)/dsqrt(psi1_sum(j)) c enddo c endif c enddo dE_old = dE_start RETURN END C ======================================================================= C SUBROUTINE FOR THE NUMERICAL INTEGRATION OF THE SCHRODINGER EQUATION C USING THE NUMEROV OR COWLING METHOD C ======================================================================= SUBROUTINE RNUMEROV(V,MD,E,Z_GRID,I_MAX,F,PSI) IMPLICIT REAL*8(A-H, O-Z) PARAMETER(npts = 1000) REAL*8 q, h, hb, k_b, eps_o, mo DATA q, h, hb, k_b, eps_o, mo 1 / 1.602D-19, 6.625D-34, 1.0544D-34, 2 1.38D-23, 8.86D-12, 9.1D-31/ REAL*8 V(0:npts) REAL*8 psi(0:npts) REAL*8 z_grid(0:npts) REAL*8 rk2(0:npts) REAL*8 norm, md C ----------------------------------------------------------------------- c calculate the k_vector and the position of the turning point c (it gives the last turning point) C ----------------------------------------------------------------------- i_match = 0 rk2(0) = 2.*md/hb*(q/hb)*(e-v(i_match)) rk = rk2(0) do i=1,i_max rk2(i)=2.*md/hb*(q/hb)*(e-v(i)) if((rk2(i)*rk.lt.0.).and.(rk.gt.0.))i_match=i rk=rk2(i) c print*,e,i,rk2(i),V(i) enddo C ----------------------------------------------------------------------- c start to generate the solution from the left C ----------------------------------------------------------------------- psimax = 0. psi(0) = 0. psi(1) = 1e-9 do i=1,i_match term1 = (z_grid(i-1)+z_grid(i))/z_grid(i-1)*psi(i) term2 = z_grid(i)/z_grid(i-1)*psi(i-1) term3 = 0.5*z_grid(i)*(z_grid(i-1)+z_grid(i))* 1 rk2(i)*psi(i) psi(i+1) = term1 - term2 - term3 if(dabs(psi(i+1)).ge.(1.D+20))then do j=1,i+1 psi(j)=psi(j)*1.D-20 enddo endif enddo PMMTCH = psi(i_match) C ----------------------------------------------------------------------- c start to generate the solution from the right C ----------------------------------------------------------------------- psi(i_max) = 0. psi(i_max-1) = 1D-9 do i=i_max-1,i_match+1,-1 term1 = (z_grid(i-1)+z_grid(i))/z_grid(i)*psi(i) term2 = z_grid(i-1)/z_grid(i)*psi(i+1) term3 = 0.5*z_grid(i-1)*(z_grid(i-1)+z_grid(i))* 1 rk2(i)*psi(i) psi(i-1) = term1 - term2 - term3 c if psi grows too large rescale all previous points if(dabs(psi(i-1)).gt.(1D+20))then do j=i_max-1,i-1,-1 psi(j)=psi(j)*1D-20 enddo endif enddo c calculation of the wavefunction that starts from right at the point c i_match-1 i = i_match term1 = (z_grid(i-1)+z_grid(i))/z_grid(i)*psi(i) term2 = z_grid(i-1)/z_grid(i)*psi(i+1) term3 = 0.5*z_grid(i-1)*(z_grid(i-1)+z_grid(i))* 1 rk2(i)*psi(i) PPMTCH = term1 - term2 - term3 C ----------------------------------------------------------------------- c renormalize the wavefunction C ----------------------------------------------------------------------- NORM = PMMTCH/psi(i_match) do i=i_match,i_max psi(i) = psi(i)*NORM enddo PPMTCH = PPMTCH*NORM psimax = dabs(psi(0)) do i=1,i_max if(dabs(psi(i)).gt.psimax)psimax=dabs(psi(i)) enddo C ----------------------------------------------------------------------- c calculate the error C ----------------------------------------------------------------------- if(psimax.ne.0)then f=(psi(i_match-1)-PPMTCH)/psimax endif RETURN END