Quantum Mechanics - The "Schroedinger" Quickbasic Program

"Schroedinger" - a Quickbasic Program to Solve Schroedinger's equation in One Dimension

Comments etc to : jo@nat.vu.nl
REM program to solve Schroedingers equation in one dimension
REM for a wavepacket incident on a localised potential feature,
REM based on Goldberg et al's paper in Am J Phys 35,177 ( 1966)
REM Written by Andrew Edgar, Physics Department, Victoria University
REM Wellington, New Zealand. e-mail Andy.Edgar@VUW.ac.nz
REM Written for the Mac using (compiled) Quickbasic.
REM
REM The potential profile is generated within the program so
REM it needs to be recompiled before it is run.
REM The internal generation of the potential profile
REM is done simply for convenience when profiles of analytic shape
REM are to be used. The program could be simply
REM modified to read in an input file specifying the potential profile
REM if required.
REM
REM The running time on a Mac IIVx with co-processor is about 50 minutes.
REM for 2000 spatial steps.
REM Be warned that the running time on other machines scales as:-
REM LCIII with coprocessor x1
REM Classic with 68030 only x6
REM Macplus x18
REM
REM ******************************************************
INPUT "# spatial steps", nx
REM  Usually use 2000  to ensure accuracy and preserve detail

REM specify the name of the output file for storing the probability amplitudes
REM
INPUT "FILENAME"; NF$
OPEN NF$ FOR OUTPUT AS #1

REM
REM starting position within range of -.5 to +.5
REM
x0=.25

REM
REM parameters defined by Goldberg et al
REM


sigma=.05 :  REM width of wavepacket
eps=1/nx: REM the spatial mesh width
eps2=eps^2
k0=3.1416/(20*eps) : REM mean k of incident wavepacket
PRINT "mean wavevector of incident packet fixed at ";k0
lambda=1: REM reflects ratio 2* eps^2/delta; delta = timestep
REM which is thus fixed by lambda
PRINT "basic step size in time="; 2*eps^2

REM wavefunction storage options:-
INPUT "Store wavefunctions only at every  nth spatial step where n=?";nmesh
stepsize = INT(.1+nmesh)
REM Wavefunction probability is calculated for all nx steps in x for accuracy
REM but only stored for x=0, nmesh*eps, 2*nmesh*eps etc. I use nmesh=8 with nx 
= 2000.
REM This is because I read the output file into Excel IV or V, and Excel can 
only
REM handle 256 columns. Furthermore, the files can become very large unless 
restricted.

REM
REM no. of time steps are now calculated
REM
nt=INT(nx*lambda/(8*eps*k0)) : REM Goldberg's condition (1)

REM
REM Similarly, only store every mtime 'th w/f otherwise files may become too 
big
REM
PRINT "# time steps="; nt
INPUT "Store wavefunctions only at every mth time step where m=?";mtime
timestep=INT(.1+mtime)

DIM V(4000),PSIR(2,4000),PSII(2,4000),eR(4000),eI(4000)
DIM fR(4000),fI(4000),omegaR(4000),omegaI(4000)

PSIR(1,0)=0
PSII(1,0)=0
PSIR(1,nx)=0
PSII(1,nx)=0

REM
REM Define the potential: this changes with every run. The
REM definition below locates a step in the centre of x=-.5 to .5
REM for 2000 steps.
REM
FOR I=0 TO 965:V(I)=0!:NEXT I
FOR I=976 TO 1035 :V(I)=1.5*(100*3.1416)^2:NEXT I
FOR I = 1036 TO 2000: V(I)=0:NEXT I

REM
REM compute and put out convergence parameters defined by Goldberg
REM

PRINT "final width increase in %="; 100*((1+(4*nt*eps2/(lambda*sigma^2))^2)^.25
-1)
PRINT "Criterion 3 (see Goldberg, should be<<1)=";(1.6*k0*eps)^2/12
REM check on the maximum value of the potential
Vmax=0:Vmin=0
FOR I= 1 TO nx STEP 1
IF V(I)>Vmax THEN Vmax = V(I)
IF V(I)<Vmin THEN Vmin=V(I)
NEXT I
PRINT "Criterion 4 (see Goldberg - should be <<1)= ";(Vmax*eps2/12)
delta=2*eps2/lambda
PRINT "Criterion 5 (see Goldberg - should be <<1) = 
";nt*(delta^3*((1.6*k0)^6-k0^6))/12
PRINT "Criterion 6 (see  Goldberg - should be <<1) =";1.6*k0*eps/3.1416
REM
REM here comes the starting form of the wavepacket, real and imaginary
REM

roll:
FOR I=1 TO nx
ex=EXP(-(I/nx-x0)^2/(2*sigma^2))
PSIR(1,I)=ex*COS(k0*I/nx)
PSII(1,I)=ex*SIN(k0*I/nx)
NEXT I

REM
REM store position
REM
PRINT #1," ";",";
FOR I=0 TO nx STEP stepsize
PRINT#1,USING "##.###";-.5+I/nx;
PRINT #1,",";
NEXT I
PRINT#1,

REM
REM store potential in output file. Note values are separated by commas.
REM
PRINT#1,"potential";",";
FOR I=0 TO nx STEP stepsize

PRINT#1, USING "##.###^^^^ ";V(I);
PRINT #1,",";
NEXT I
PRINT#1,

REM
REM store psi psi* for first frame
REM note max value for later scale choice
REM

probmax=0

PRINT#1,"time=0";",";
FOR I=0 TO nx STEP stepsize
prob=(PSIR(1,I)^2+PSII(1,I)^2)
PRINT#1, USING "##.#######"; prob;
PRINT  #1,",";
IF probmax<prob THEN probmax=prob
NEXT I
PRINT#1,

REM
REM function e defined by Goldberg
REM

eR(1)=2+eps2*V(1)
eI(1)=-lambda
FOR I=2 TO nx
e2=eR(I-1)^2+eI(I-1)^2
eR(I)=2+eps2*V(I)-eR(I-1)/e2
eI(I)=-lambda+eI(I-1)/e2
NEXT I

REM
REM main loop
REM omega, f,psi with real and imaginary parts as defined by Goldberg et al;
REM

FOR j=1 TO  nt STEP timestep
PRINT "CALCULATING - STEP ";timestep+(j-1) : REM for the impatient, announce 
how far you've got
FOR k=1 TO timestep
FOR I= 1 TO nx-1
omegaR(I)=-PSIR(1,I+1)-PSIR(1,I-1)+(2+eps2*V(I))*PSIR(1,I)-lambda*PSII(1,I)
omegaI(I)=-PSII(1,I+1)-PSII(1,I-1) + lambda*PSIR(1,I) +(2+eps2*V(I))*PSII(1,I)
NEXT I
fR(1)=omegaR(1)
fI(1)=omegaI(1)
FOR I=2 TO nx-1
e2=eR(I-1)^2+eI(I-1)^2
fR(I)=omegaR(I)+(fR(I-1)*eR(I-1)+fI(I-1)*eI(I-1))/e2
fI(I)=omegaI(I)+(fI(I-1)*eR(I-1)-fR(I-1)*eI(I-1))/e2
NEXT I
PSIR(2,nx-1)=(-fR(nx-1)*eR(nx-1)-fI(nx-1)*eI(nx-1))/e2
PSII(2,nx-1)=(+fR(nx-1)*eI(nx-1)-fI(nx-1)*eR(nx-1))/e2
FOR I=nx-2 TO 1 STEP -1
e2=eR(I)^2+eI(I)^2
PSIR(2,I)=((PSIR(2,I+1)-fR(I))*eR(I)+(PSII(2,I+1)-fI(I))*eI(I))/e2
PSII(2,I)=((PSII(2,I+1)-fI(I))*eR(I)-(PSIR(2,I+1)-fR(I))*eI(I))/e2
NEXT I
FOR I=0 TO nx
PSIR(1,I)=PSIR(2,I)
PSII(1,I)=PSII(2,I)
NEXT I
NEXT k

REM
REM save psi psi*
REM

PRINT #1,"time=";
PRINT#1,USING "#.###^^^^";2*eps^2*(timestep+(j-1));
PRINT#1,",";

FOR I=0 TO nx STEP stepsize
prob=(PSIR(1,I)^2+PSII(1,I)^2)
PRINT#1, USING "##.########"; prob;
PRINT  #1,",";
IF probmax<prob THEN probmax=prob
NEXT I
PRINT#1,
NEXT j
PRINT #1, "max pot.="; Vmax
PRINT #1, "min pot.="; Vmin
PRINT #1, "max prob.="; probmax
PRINT #1,"nx=";nx
PRINT #1,"no time steps=";nt
PRINT #1, "sigma=";sigma
PRINT #1, "epsilon=";eps
PRINT #1,"k0=";k0
PRINT #1,"x step size for save=";stepsize
PRINT #1,"t step size for save=";timestep
CLOSE #1
STOP
END