C PROGRAM SCHR1D C The following computer code provides the solution of the 1-D C time-dependent Schroedinger equation using the Crank-Nicholson C semi--implicit approach. The initial condition is a gaussian C a gaussian wavepacket with a specified kinetic energy C corresponding to the motion of the packet centroid. At the C boundaries the wavefunction is assumed to be zero, which is C equivalent to perfectly reflecting boundary conditions. C The code reads three input files which must be provided by C the user: C file INPUT1.D, which includes a number of simulation parameters C file EFMASS.D, which contains the relative effective mass C defined on the mesh intervals C (so that EFMASS(i) = $m^*(i-1/2)$) C file VOLT.D, with the potential energy distribution in Joule C units, defined on the mesh points. C Basic version in complex arithmetic of the 1-D time-dependent C Schroedinger equation, with space-dependent effective mass. C Crank-Nicholson finite difference approach. Subroutine C TRIDIAG solves the resulting tridiagonal matrix. C Initial condition is a Gaussian wavepacket moving from C left to right. C No scaling used - M.K.S. units. Energy input in eV units. C---------------------------------------------------------------------- C Definitions REAL K0 COMPLEX PSI,C2,UIM,BC1,BC2,A,B,C,RHS,W,G C---------------------------------------------------------------------- C COMMMON blocks C V = potential energy profile C PSI = wavefunction amplitude (complex) C PACK = wavefunction squared (real) C EEE = variable for generation of gaussian C EFMASS = space dependent relative effective mass C EMINV1 = inverse of relative effective mass C EMINV2 = sum of adjacent vales of effective mass inverse C A,B,C,RHS = coefficients for matrix solution C W,G = working arrays for matrix solution C N = number of mesh points C COMMON /VAR1/ V(20000),PSI(20000),PACK(20000) COMMON /VAR2/ EEE(20000) COMMON /VAR3/ EFMASS(20000),EMINV1(20000),EMINV2(20000) COMMON /MAT1/ A(20000),B(20000),C(20000),RHS(20000) COMMON /MAT2/ W(20000),G(20000) COMMON /MESH/ N C---------------------------------------------------------------------- C Open input/output files OPEN(unit=8,file='INPUT1.D') OPEN(unit=9,file='EFMASS.D') OPEN(unit=10,file='VOLT.D') OPEN(unit=20,file='OUT1.DAT') OPEN(unit=30,file='OUT2.DAT') C---------------------------------------------------------------------- C Physical parameters C HBAR=Planck's constant/(2*pi) C ELMASS = electron mass in vacuum C---------------------------------------------------------------------- HBAR = 1.0546E-34 ELMASS = 9.1E-31 C---------------------------------------------------------------------- C Imaginary unit UIM=(0.,1.) C---------------------------------------------------------------------- C INPUT - Read discretization parameters and energy C NITER = total number of time iterations C DX = mesh interval C DT = timestep C SIGFACT = spread of initial wavepacket in units of DX C VVT = energy of particle [eV] C N = number of internal mesh points (excluding boundaries) C ICENT = initial mesh point location of wavepacket centroid C---------------------------------------------------------------------- READ(8,*) NITER READ(8,*) DX READ(8,*) DT READ(8,*) SIGFACT READ(8,*) VVT READ(8,*) N READ(8,*) ICENT C---------------------------------------------------------------------- C INPUT - Effective mass distribution DO 10 I=1,N READ(9,*) EFMASS(I) 10 CONTINUE C---------------------------------------------------------------------- C INPUT - Potential distribution DO 20 I=1,N READ(10,*) V(I) 20 CONTINUE C---------------------------------------------------------------------- C Initialize time clock TIME = 0.0 C---------------------------------------------------------------------- C Simple zero boundary conditions for wavefunction BC1=(0.,0.) BC2=(0.,0.) C---------------------------------------------------------------------- C Form initial gaussian wavepacket C K0 = momentum of wavepacket centroid, obtained from electron C kinetic energy for parabolic band C SIGFACT = gaussian spread in units of DX C SIG = resulting gaussian spread (r.m.s.) C---------------------------------------------------------------------- K0 = SQRT(2.*EFMASS(ICENT)*ELMASS/HBAR*VVT/HBAR) SIG = SIGFACT*DX C---------------------------------------------------------------------- DO 90 I = 1,N EEE(I) = EXP(-((I)*DX-ICENT*DX)**2/2./SIG**2) PSI(I) = COS(K0*(I)*DX)*EEE(I)*10+ * SIN(K0*(I)*DX)*EEE(I)*10*UIM 90 CONTINUE C---------------------------------------------------------------------- C SUM0 = Simple Integral over the domain of wavepacket probability C---------------------------------------------------------------------- SUMO=0. DO 100 I = 1,N SUMO=SUMO+ABS(PSI(I))**2 100 CONTINUE SUMO = SUMO*DX C---------------------------------------------------------------------- C Normalize wavepacket DO 400 I = 1,N PSI(I) = PSI(I) / SQRT(SUMO) 400 CONTINUE DO 500 I = 1,N PACK(I) = ABS(PSI(I))**2 500 CONTINUE WRITE (20,*) TIME,SUMO DO 550 I=1,N WRITE (20,*) I,PACK(I) 550 CONTINUE C---------------------------------------------------------------------- C Set up coefficients for matrix calculation C---------------------------------------------------------------------- C1 = DX**2/HBAR*2.*ELMASS/HBAR C2 = (2.*DX**2/DT*2.*ELMASS/HBAR)*UIM C---------------------------------------------------------------------- DO 600 I=1,N-1 EMINV1(I) = 1./EFMASS(I) 600 CONTINUE DO 650 I=1,N-1 EMINV2(I) = EMINV1(I) + EMINV1(I+1) 650 CONTINUE C---------------------------------------------------------------------- C DO 1000 - OUTER TIME LOOP C---------------------------------------------------------------------- DO 1000 ITER=1,NITER C---------------------------------------------------------------------- C Evaluation of coefficients of the tridiagonal matrix C B(i) = Main diagonal C A(i) = Lower diagonal C C(i) = Upper diagonal C---------------------------------------------------------------------- DO 200 I=1,N B(I)=-EMINV2(I)+C2-C1*V(I) C(I)= EMINV1(I+1) A(I)= EMINV1(I) 200 CONTINUE A(1)=0. C(N)=0. C---------------------------------------------------------------------- C Evaluation of the right hand side C---------------------------------------------------------------------- DO 300 I=1,N RHS(I)=-PSI(I-1)*EMINV1(I)+(C2+C1*V(I)+EMINV2(I))*PSI(I) # -PSI(I+1)*EMINV1(I+1) 300 CONTINUE RHS(1)=RHS(1)-BC1 RHS(N)=RHS(N)-BC2 C---------------------------------------------------------------------- CALL TRIDIAG C---------------------------------------------------------------------- C PACK = Wavefunction probability (squared modulus) C---------------------------------------------------------------------- DO 5000 I = 1,N PACK(I) = ABS(PSI(I))**2 5000 CONTINUE C---------------------------------------------------------------------- C Check integral of wavefunction C---------------------------------------------------------------------- SUM = 0. DO 6000 I = 1,N SUM = SUM+PACK(I)*DX 6000 CONTINUE C---------------------------------------------------------------------- C Update time clock TIME = TIME+DT C---------------------------------------------------------------------- 1000 CONTINUE C---------------------------------------------------------------------- C END OF TIME LOOP C---------------------------------------------------------------------- WRITE (30,*) TIME,SUM DO 3000 I=1,N WRITE (30,*) I, PACK(I) 3000 CONTINUE C---------------------------------------------------------------------- C Close input/output files C---------------------------------------------------------------------- CLOSE(8) CLOSE(9) CLOSE(10) CLOSE(20) CLOSE(30) C---------------------------------------------------------------------- STOP END C---------------------------------------------------------------------- SUBROUTINE TRIDIAG COMMON /VAR1/ V(20000),PSI(20000),PACK(20000) COMMON /MAT1/ A(20000),B(20000),C(20000),RHS(20000) COMMON /MAT2/ W(20000),G(20000) COMMON /MESH/ N COMPLEX PSI,A,B,C,RHS,W,G C---------------------------------------------------------------------- C Gaussian Elimination C---------------------------------------------------------------------- W(1)=C(1)/B(1) DO 10 I=2,(N-1) W(I)=C(I)/(B(I)-A(I)*W(I-1)) 10 CONTINUE G(1)=RHS(1)/B(1) DO 20 I=2,N G(I)=(RHS(I)-A(I)*G(I-1))/(B(I)-A(I)*W(I-1)) 20 CONTINUE C---------------------------------------------------------------------- C Backsubstitution C---------------------------------------------------------------------- PSI(N)=G(N) DO 30 I=1,(N-1) K=N-I PSI(K)=G(K)-W(K)*PSI(K+1) 30 CONTINUE C---------------------------------------------------------------------- RETURN END