program tdse implicit none ! Motion of wavepacket incident on a potential integer n_total parameter (n_total=2100) real x0,k0,width,V0,a,xmin,xmax,dx,dx2,t,dt,V1 real Re(n_total),Im(n_total),Imold(n_total) real Re0(n_total),Im0(n_total),PP0(n_total) integer i, n , s ,s0 , s1 , p1 , p2 , p3 , p4 integer pgopen p1=1 p2=2 p3=3 p4=4 call parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt) c if (n.gt.n_total) stop ('Not enough Grids in n_total') ! Not for g77 if (n.gt.n_total) then write(*,*) 'Not enough Grids in n_total' end if call initial_packet(Re,Im,x0,k0,width,xmin,n,dx,dt) if (pgopen('/xwin').le.0) stop call set_up_windows(xmin,xmax,V0,p1,p2,p3,p4) call draw_potential_2(V0,a,xmin,xmax,p4) c Re0=Re ! Array operation is dangerous to use in some versions. c Im0=Im c do i=1,n c Re0(i)=Re(i) c Im0(i)=Im(i) c end do s0=0 s1=0 do while (s0==0) s1=0 c do while(s1==0) c open (UNIT=10,FILE='control',STATUS='OLD',IOSTAT=s1) c if (s1>0) exit c read(10,*) s1 c close(10) c end do c call delay_loop(10000000) call delay_loop(100000) call evolve(Re,Im,Imold,t,V0,a,dx,dx2,dt,xmin,n) call plots(Re,Re0,Im,Im0,Imold,PP0,t,xmin,n,dx,p1,p2,p3,p4) end do stop end program tdse subroutine parameters(x0,k0,width,V0,a,xmin,xmax,n,dx,dx2,dt) implicit none real x0,k0,width,V0,a,dx,dx2,dt real xmin , xmax integer n x0 = -50.0 c width =4.0 width=1.0 k0 = 2.0 xmax = 70.0 xmin = -1.0*xmax V0 = 2.0 a = 0.5 dx = 0.40 dx2 = dx*dx n = int((xmax-xmin)/dx) + 1 dt = 0.1 c dt = 0.01 c Check V0 based on Visscher's criteria for stable algorithm. if ( (V0.lt.(-2.0/dt)) .OR. (V0.gt.2.0*(1/dt - 1/dx2))) then write(*,*) 'Warning, V0 too big or too small.' end if return end ! initial Gaussian wavepacket subroutine initial_packet(Re,Im,x0,k0,width,xmin,n,dx,dt) implicit none real x0,k0,width,xmin,dx,dt,A,b,arg,xmax real delta2, x, t, e real Re(*), Im(*) real A0,A1,A2 real width_1 , width_2 real pi integer i, n A1=2.0 A2=2.0 width_1=4.0 width_2=4.0 pi = 3.1415926 width=width_1 delta2 = width *width A = (2*pi*delta2)**(-0.25) b = 0.5*k0*dt A0=A A=A0*A1 do i = 1 , n x = xmin + (i-1)*dx arg = 0.25 * ((x - x0)**2)/delta2 e = exp(-arg) ! value at t = 0 Re(i) = A*cos(k0*(x - x0))*e !value at t = 0 ! value at t = dt/2 Im(i) = A*sin(k0*(x - x0 - 0.5*b))*e !value at t = dt/2 end do goto 222 ! Skip second wave packet generation c Second Wave Packet A=A0*A2 width=width_2 delta2 = width *width x0= -1*x0 k0= -1*0.5*k0 b = 0.5*k0*dt do i = 1 , n x = xmin + (i-1)*dx arg = 0.25 * ((x - x0)**2)/delta2 e = exp(-arg) ! value at t = 0 Re(i) = Re(i) + A*cos(k0*(x - x0))*e !value at t = 0 ! value at t = dt/2 Im(i) = Im(i) + A*sin(k0*(x - x0 - 0.5*b))*e !value at t = dt/2 end do 222 continue return end subroutine set_up_windows(xmin,xmax,V0,p1,p2,p3,p4) implicit none integer pgopen, p1, p2, p3, p4 real V0, xmin, xmax ,V1 real x_set(2),y_set(2) p1=1 p2=2 p3=3 p4=4 x_set(1)=xmin y_set(1)=0.0 x_set(2)=xmax y_set(2)=0.0 ! if (pgopen('/xwin').le.0) stop call pgpap (4.0,2.0) call pgsubp (1,4) c call pgenv (xmin,xmax,V1,V0,0,0) c call pgpanl(1,1) call pgenv (xmin,xmax,-1.0,4.0,0,0) call pgenv (xmin,xmax,-1.0,1.0,0,0) c call pgline(2,x_set,y_set) c call pgpanl(1,2) call pgenv (xmin,xmax,-1.5,1.0,0,0) c call pgline(2,x_set,y_set) c call pgpanl(1,3) call pgenv (xmin,xmax,-1.0,1.0,0,0) c call pgline(2,x_set,y_set) V1=-1.0*V0 c call pgpanl(1,4) c call pgenv (xmin,xmax,V1,V0,0,0) c call pgsvp(xmin,xmax,-1.0,1.0) return end subroutine draw_potential_2 (V0,a,xmin,xmax,p4) implicit none real V0, a, xmin, xmax, delta_x, V integer i, p4, n_grids parameter (n_grids=1000) real xx(n_grids), yy(n_grids), ymin, ymax, y_range external V delta_x = (xmax-xmin)/(n_grids-1) do i=1, n_grids xx(i) = xmin + (i-1)*delta_x end do do i=1, n_grids yy(i) = V(xx(i),V0,a) end do ymin = yy(1) ymax = yy(1) do i=2, n_grids if(yy(i).lt.ymin) ymin = yy(i) if(yy(i).gt.ymax) ymax = yy(i) end do y_range = ymax - ymin y_range = 1.2*y_range do i=1, n_grids yy(i) = yy(i)/y_range end do call pgpanl(1,p4) call pgsci(2) call pgline(n_grids,xx,yy) call pgsci(1) end subroutine plots(Re,Re0,Im,Im0,Imold,PP0,t,xmin,n,dx,p1, &p2,p3,p4) implicit none real Re(*),Im(*),Imold(*) real Im0(*),Re0(*),PP0(*) real t,xmin,dx,xmax integer n_plot_max parameter (n_plot_max=1100) real Psum , P , x , xx(n_plot_max) , PP(n_plot_max) real x_plots(2) , y_plots(2) integer p1,p2,p3,p4,i,n if (n_plot_max.lt.n) then write(*,*) 'Not enough grid for ploting, try increase n_plot.' endif p1=1 p2=2 p3=3 p4=4 xmax=-1.0*xmin x_plots(1)=xmin y_plots(1)=0.0 x_plots(2)=-1*xmin y_plots(2)=0.0 call pgpanl(1,p2) call pgsci(3) call pgline(2,x_plots,y_plots) call pgsci(1) do i = 1, n x = xmin + (i-1)*dx xx(i)=x end do call pgbbuf call pgsci(0) call pgline(n,xx,Re0) call pgsci(1) c call pgenv(xmin,xmax,-1.0,1.0,0,-1) c call pgsvp(xmin,xmax,-1.0,1.0) call pgline(n,xx,Re) call pgebuf do i =1, n Re0(i)=Re(i) end do call pgpanl(1,p3) call pgsci(3) call pgline(2,x_plots,y_plots) call pgsci(1) call pgbbuf call pgsci(0) call pgline(n,xx,Im0) call pgsci(1) call pgline(n,xx,Im) call pgebuf do i =1 , n Im0(i)=Im(i) end do call pgpanl(1,p1) call pgsci(3) call pgline(2,x_plots,y_plots) call pgsci(1) Psum = 0.0 do i = 1, n x = xmin + (i-1)*dx P = Re(i)*Re(i) + Im(i)*Imold(i) PP(i)=P Psum = Psum + P*dx end do call pgbbuf call pgsci(0) call pgline(n,xx,PP0) ! call pgeras ! call pgpanl(1,1) ! call pgenv (xmin,xmax,-1.0,4.0,0,0) call pgsci(1) ! call pgenv (xmin,xmax,-1.0,4.0,0,0) call pgline(n,xx,PP) call pgebuf do i =1 , n PP0(i)=PP(i) end do return end subroutine evolve(Re,Im,Imold,t,V0,a,dx,dx2,dt,xmin,n) implicit none c real Re(*), Im(*), Imold(*) real Re(1100), Im(1100), Imold(1100) real t,V0,a,dx,dx2,dt,xmin,HRe,HIm,x,V,xmax external V integer n, i c xmax=-1.0*xmin do i = 1 , n x=xmin+(i-1)*dx c real part defined at multiples of dt HIm= V(x,V0,a)*Im(i)-0.5*(Im(i+1)-2*Im(i)+Im(i-1))/dx2 Re(i) = Re(i)+HIm*dt end do do i = 1 , n x=xmin+(i-1)*dx c dt/2 earlier than real part Imold(i)=Im(i) HRe=V(x,V0,a)*Re(i)-0.5*(Re(i+1)-2*Re(i)+Re(i-1))/dx2 c dt/2 later than real part Im(i) = Im(i) - HRe*dt end do t = t + dt return end function V(x,V0,a) implicit none real x, V0, a real V c Potential Step c if (x .gt. a) then c V=V0 c else c V=0 c end if c Potential Barrier V = 0.0 if ( (x.gt.(-1*a)) .AND. (x.lt.a) ) V = V0 c Infinite potential well c V = 100000.0 c a = 50.0 c if ( (x.gt.(-1*a)) .AND. (x.lt.a) ) V = 0.0 end subroutine delay_loop (m) real a integer i, m do i=1, m a=3.14159 a= sqrt (a) end do return end