C--------------------------------------------------------------- C C C PHYSICS ON PARALLEL COMPUTERS - PROJECT C C "Scattering of Quantum Wave Packets" C C Peter Hammerstein C C Abstract: C C The Time Dependent Schrodinger Equation is solved numerically in order C to investigate the scattering of Quantum Wave Packets from various time C independent potentials. In general such problems are not solvable analytically. C---------------------------------------------------------------- implicit none include 'pict_subs.h' real, parameter:: pi = 3.1415927 real, parameter:: mass = 1.0 real, parameter:: h_bar= 1.0 integer, parameter:: bits = 8 integer, parameter:: length=256 ! latice length integer i, ii, j, step character decide real Xo, Yo, Ko_x, Ko_y ! initial wave values real sigma,t, t_step, x_unit ! values for wave packet real kin_energy,pot_energy, psisquare ! functions real expected_x, expected_y ! functions complex alpha, beta logical A, B, C, D, E, F, G, H complex, dimension(length,length):: psi, pot, e_pot ! wave & pot. complex, dimension(length,length):: betaeven, betaodd complex, dimension(length,length):: visual ! display array !_for potential_____________ character choice real pot_value integer help1, help2, help3, help4, help5 real help6 logical, dimension(length, length):: zyl !_for plot__________________ integer rnx, rny, nxp, nyp, display C_____ALL ARRAYS ARE LOCATED ON PARALLEL PROCESSORS_________________ CMF$ layout PSI(:news, :news) CMF$ layout POT(:news, :news) CMF$ layout betaeven(:news, :news) CMF$ layout betaodd(:news, :news) CMF$ layout e_pot(:news, :news) CMF$ layout visual(:news, :news) C___INITIALISATION: PARAMETERS FOR INTEGRATION, POTENTIAL AND WAVE____ 8 H = .false. A = .false. B = .false. C = .false. D = .false. E = .false. F = .false. G = .false. write(*,*) ' ' write(*,*)'____________________________________________________' write(*,*)'TIME STEP dt: ' read(*,*) t_step write(*,*)'SPACE UNIT dx: ' read(*,*) x_unit write(*,*)'A QUANTUM WAVE PACKET WILL BE SCATTERED AT A' write(*,*)'POTENTIAL OF YOUR CHOISE. ENTER YOUR LETTER:' 16 write(*,*)' ' write(*,*)' YLINDRICAL POTENTIAL' write(*,*)' ARRIER ' write(*,*)' SOT BARRIER ' write(*,*)' LOT POTENTIAL ' write(*,*)' IRCLE BARRIER ' write(*,*)' ARMONIC POTENTIAL' write(*,*)' O POTENTIAL' read(*,*) choice if (choice.ne.'N'.and.choice.ne.'S'.and.choice.ne.'B' & .and.choice.ne.'Z'.and.choice.ne.'H'.and.choice.ne.'F' & .and.choice.ne.'C') then write(*,*)'!!! C A P I T A L L E T T E R S !!!' goto 16 end if !____________________________________________________________________ !__________CONFIGURATION OF THE POTENTIAL____________________________ pot = (0.0,0.0) if (choice.ne.'N') then write(*,*)'ENTER EXTREME VALUE FOR THE POTENTIAL (! sign !):' read(*,*) pot_value else goto 32 end if C______HARMONIC POTENTIAL_____________________________ if (choice.eq.'H') then write(*,*)'ENTER POTENTIAL CENTRE:' read(*,*) help1, help2 ! centre_x, centre_y write(*,*)'ENTER RADIUS:' read(*,*) help6 ! radius forall(i=1:length,j=1:length) & pot(i,j) = ((sqrt(abs(pot_value))/help6 * (i-help1))**2. + & (sqrt(abs(pot_value))/help6 * (j-help2))**2. - & abs(pot_value)) * pot_value/abs(pot_value)*(-1) if (pot_value.gt.0) then where (real(pot).lt.0) pot = 0.0 end if if (pot_value.lt.0) then where(real(pot).gt.0) pot = 0.0 end if where (real(pot).ne.0.0) pot = pot - pot_value where (real(pot).ne.0.0) zyl = .true. elsewhere zyl = .false. end where end if C_______SHARP ZYLINDRICAL POTENTIAL_______________________ if (choice.eq.'Z') then write(*,*)'ENTER POTENTIAL CENTRE :' read(*,*) help1, help2 ! centre_x, centre_y write(*,*)'ENTER RADIUS:' read(*,*) help3 ! radius forall(i=1:length,j=1:length) & zyl(i,j) = (((i-help1)**2. + (j-help2)**2.).le.help3**2.0) where (zyl) pot = pot_value end if C_______BARRIER PARALLEL Y-AXIS_________________________ if (choice.eq.'B') then 24 write(*,*)'ENTER X_START AND X_END VALUES FOR BARRIER :' read(*,*) help1, help2 ! barrier_start, barrier_end if ((help1.ge.help2).or.(max(help1,help2).gt.length)) goto 24 forall(i=help1:help2) POT(i,:) = pot_value end if C______SOFT BARRIER PARALLEL Y-AXIS_____________________ if (choice.eq.'F') then 28 write(*,*)'ENTER X_START AND X_END VALUES FOR BARRIER :' write(*,*)' !!!! (X_END - X_START) > 3 !!!!' read(*,*) help1, help2 ! barrier_start, barrier_end if ((help1.ge.help2).or.(max(help1,help2).gt.length)) goto 28 forall(i=help1:help2) POT(i,:) = pot_value forall(i=help1-10:help1+10) & pot(i,:) = (1.-(help1+10-i)*0.05)*pot_value forall(i=help2-10:help2+10) & pot(i,:) = (1.-(i-(help2-10))*0.05)*pot_value end if C_______SLOTS_______________________________________ if (choice.eq.'S') then write(*,*)'ENTER NUMBER OF SLOTS :' read(*,*) help1 ! number of slots write(*,*)'ENTER WIDTH OF SLOT (ALL ARE THE SAME)' read(*,*) help2 write(*,*)'ENTER X_START AND X_END OF ''WALL'':' read(*,*) help3, help4 forall(i=help3:help4,j=1:length) pot(i,j) = pot_value do ii=1,help1 forall(i=help3:help4,j=ii*((length-help1*help2)/(help1+1)) & + (ii-1)*help2 + 1 : ii*(((length-help1*help2)/ & (help1+1))+help2)) pot(i,j) = 0.0 end do end if C_______CIRCLE BARRIER__________________________________ if (choice.eq.'C') then write(*,*)'ENTER CENTRE OF CIRCLE' read(*,*) help1, help2 write(*,*)'ENTER INNER AND OUTER RADIUS' read(*,*) help3, help4 forall(i=1:length,j=1:length) & zyl(i,j) = (((i-help1)**2. + (j-help2)**2.).le.help4**2.0) where (zyl) pot = pot_value forall(i=1:length,j=1:length) & zyl(i,j) = (((i-help1)**2. + (j-help2)**2.).le.help3**2.0) where (zyl) pot = 0.0 end if C______________________________________________________________________ C______PLOT POTENTIAL DATA TO FILE FOR EXTERNAL VISUALISATION_______ write(*,*)'DO YOU WHISH EXTERNAL PLOT OF POTENTIAL ?' read(*,*) decide if (decide.ne.'n') then H = .true. open(unit=10,file='pot.dat',status='unknown') do i=1,length write(10,*) i, real(pot(i,128)) end do end if c______PLOT SOME POTENTIAL DATA FOR SCREEN CONTROL___________________ write(*,*)' X potential ' do i=1,length,10 write(*,*) i, real(pot(i,128)) end do write(*,*)'IS EVERYTHING SATISFYING YOU , OR ???' read(*,*) decide if (decide.eq.'n') goto 72 C______SHOW DISPLAY WINDOW___________________________________________ 32 call init_display(length, length, bits) C_______VALUES FOR INITIAL WAVE______________________________________ write(*,*)' ' write(*,*) 'YOU WILL USE A GAUSSIAN WAVE PACKET !!!' write(*,*) 'ENTER STANDARD DEVIATION (SIGMA):' read (*,*) sigma sigma = x_unit * sigma write(*,*) 'ENTER INITIAL CENTRE OF WAVE PACKET (Xo, Yo)' read (*,*) Xo, Yo Xo = x_unit * Xo Yo = x_unit * Yo write(*,*) 'ENTER START WAVE VECTOR(K0_x,K0_y)' read (*,*) Ko_X, Ko_Y C_____SET INITIAL WAVE______________________________________________ forall(i=1:length,j=1:length) psi(i,j) = & 1.0 / (sqrt( 2.0 * pi )* sigma) & * exp( (0.,1.) * x_unit * (Ko_X * i + Ko_Y * j) ) & * exp( -(((x_unit * i - Xo)**2.0 & + (x_unit * j - Yo)**2.0) / (4.0*sigma**2.0))) !!! INITIAL WAVE IS SCALED WITH X_UNIT !!! C_______DISPLAY INITIAL WAVE (PSI) AND POTENTIAL ____________________ visual = psi where (pot.ne.0.0) visual = maxval(abs(psi)) call display_psi(visual/maxval(abs(psi)), length, length) write(*,*) 'psi square = ', psisquare(length,psi,x_unit) write(*,*)'IS EVERYTHING SATISFYING YOU , OR ???' read(*,*) decide if (decide.eq.'n') goto 72 C_____CALCULATE ALPHA AND BETA FOR "RICHARDSON ALGORITHMUS"__________ alpha = 0.5 * (1+exp( (0.0,-1.0) * t_step * h_bar / & (mass * x_unit**2.0)) ) beta = 0.5 * (1-exp( (0.0,-1.0) * t_step * h_bar / & (mass * x_unit**2.0)) ) C_____BETA MASKS___________________________ forall(i=1:length,j=1:length) & betaeven(i,j) = beta * (1.-(-1.)**(i+j))/2. forall(i=1:length,j=1:length) & betaodd(i,j) = beta * (1.+(-1.)**(i+j))/2. C_____HELP VALUE FOR SUBROUTINE_________ e_pot = exp( (0.0,-1.0) * t_step * pot / h_bar) step = 0 C_____SEVERAL VALUES CAN BE PRINTED TO FILES__________________________ write(*,*)' ' a = .true. write(*,*)'wish data printed to file ?' read(*,*) decide if (decide.ne.'n') then write(*,*)' expected x over step ???' read(*,*) decide if (decide.ne.'n') A = .true. write(*,*)' expected y over step ???' read(*,*) decide if (decide.ne.'n') B = .true. write(*,*)' kin. Energy over step ???' read(*,*) decide if (decide.ne.'n') C = .true. write(*,*)' pot Energy over step ???' read(*,*) decide if (decide.ne.'n') D = .true. write(*,*)' total Energy over step ???' read(*,*) decide if (decide.ne.'n') E = .true. write(*,*)' psi**2 over step ???' read(*,*) decide if (decide.ne.'n') F = .true. write(*,*)' for circle only: psi**2 inside circle over step' read(*,*) decide if (decide.ne.'n') G = .true. end if if (A) open(unit=20,file='expect_x.dat',status='unknown') if (B) open(unit=30,file='expect_y.dat',status='unknown') if (C) open(unit=40,file='ekin.dat',status='unknown') if (D) open(unit=50,file='epot.dat',status='unknown') if (E) open(unit=60,file='etot.dat',status='unknown') if (F) open(unit=70,file='abs.dat',status='unknown') if (G) open(unit=80,file='decay.dat',status='unknown') write(*,*) ' ' write(*,*)'' read(*,*) !_____________________________________________________________________ C________M A I N P A R T____________________________________________ do while(1.le.5) step = step + 1 t = t + t_step call psi_update(length,alpha,betaeven,betaodd,psi,e_pot) C______RUN MANIPULATION_________________ if (mod(step,500).eq.0) then 64 write(*,*)'how with potential ? ; continue the computation & , es or o ???' read(*,*) decide if (decide.eq.'s') then visual = psi where(pot.ne.0.0) visual = maxval(abs(psi)) call display_psi(visual/maxval(abs(psi)), length, length) write(*,*)'how without potential ? ; continue the & computation , es or o ???' read(*,*) decide if (decide.eq.'s') then call display_psi(psi/maxval(abs(psi)), length, length) goto 64 end if end if if (decide.eq.'n') goto 72 end if C__________DISPLAY PSI_____________________ if (mod(step,50).eq.0) then call display_psi(psi/maxval(abs(psi)), length, length) end if C_________SCREEN AND FILE OUTPUT OF DIFFERENT VALUES_____ if (mod(step,10).eq.0) then write(*,*) '*' write(*,*)'step = ',step write(*,*)'psisquare = ',psisquare(length, psi,x_unit) write(*,*)'E_kin = ', kin_energy(length,psi, pot) write(*,*)'E_pot = ',pot_energy(length, x_unit, psi, pot), write(*,*)'E_total = ', kin_energy(length,psi, pot) + & pot_energy(length,x_unit,psi,pot) write(*,*)'expected x = ',expected_x(psi, x_unit, length) write(*,*)'expected y = ',expected_y(psi, x_unit, length) if (A) write(20,*) step, expected_x(psi,x_unit,length) if (B) write(30,*) step, expected_y(psi,x_unit,length) if (C) write(40,*) step, kin_energy(length,psi,pot) if (D) write(50,*) step, pot_energy(length,x_unit,psi,pot) if (E) write(60,*) step, kin_energy(length,psi, pot) + & pot_energy(length,x_unit,psi,pot) if (F) write(70,*) step, psisquare(length, psi, pot) if (G) write(80,*) step, sum((abs(psi)**2),mask=zyl)*x_unit**2 end if end do ! END MAIN PART 72 write(*,*)'!!! new try, new luck !!!' if (H) close(10) if (A) close(20) if (B) close(30) if (C) close(40) if (D) close(50) if (E) close(60) if (F) close(70) if (G) close(80) C______LAST SCREEN OUTPUT______________________________________________ write(*,*)' ' write(*,*)'dt= ',t_step,' x_unit= ',x_unit,' sigma= ',sigma write(*,*)'potential type was : ',choice write(*,*)' ' 80 write(*,*)'gain or eave ?' read(*,*) decide if (decide.eq.'a') goto 8 if (decide.ne.'a'.and.decide.ne.'l') goto 80 stop end C_________S U B R O U T I N E ________________________________ C_________RICHARDSON ALGORITHM:_______________________ subroutine psi_update(length,alpha,betaeven,betaodd,psi,e_pot) implicit none integer length complex alpha complex, dimension(length,length):: psi, e_pot complex, dimension(length,length):: betaeven, betaodd cmf$ layout psi(:news, :news) cmf$ layout e_pot(:news, :news) cmf$ layout betaeven(:news,:news) cmf$ layout betaodd(:news,:news) psi = psi * e_pot psi = alpha * psi + betaeven * cshift(psi,dim=1,shift=-1) & + betaodd * cshift(psi,dim=1,shift=+1) psi = alpha * psi + betaeven * cshift(psi,dim=1,shift=+1) & + betaodd * cshift(psi,dim=1,shift=-1) psi = alpha * psi + betaeven * cshift(psi,dim=2,shift=-1) & + betaodd * cshift(psi,dim=2,shift=+1) psi = alpha * psi + betaeven * cshift(psi,dim=2,shift=+1) & + betaodd * cshift(psi,dim=2,shift=-1) return end C__________ F U N C T I O N S _____________________________ C_______NORM OF THE WAVE PACKET_________ real function psisquare(length, array, x_unit) implicit none integer length real x_unit complex, dimension(length,length):: array cmf$ layout array(:news, :news) psisquare = sum(abs(array)**2) * x_unit**2 return end C_________KINETIC ENERGY__________________ real function kin_energy(length,psi) implicit none integer length real, parameter:: h_bar = 1.0 real, parameter:: mass = 1.0 complex, dimension(length,length)::psi cmf$ layout psi(:news, :news) kin_energy = & - h_bar**2/(2*mass)*sum(conjg(psi)* & ( cshift(psi, dim=1, shift=+1) + & cshift(psi, dim=1, shift=-1) + & cshift(psi, dim=2, shift=+1) + & cshift(psi, dim=2, shift=-1) - 4*psi) ) return end C__________POTENTIAL ENERGY_________________ real function pot_energy(length,x_unit,psi, pot) implicit none real x_unit integer length complex, dimension(length,length)::psi real, dimension(length,length):: pot cmf$ layout psi(:news, :news) cmf$ layout pot(:news,:news) pot_energy = sum(abs(psi)**2 * pot ) * x_unit**2 return end C_________EXPECTED X-VALUE___________________ real function expected_x(psi,x_unit,length) implicit none integer i, j, length real x_unit complex, dimension(length,length):: psi,xpsisquare CMF$ LAYOUT psi(:news,:news) CMF$ LAYOUT xpsisquare(:news,:news) forall(i=1:length,j=1:length) & xpsisquare(i,j) = i * x_unit * abs(psi(i,j))**2 expected_x = sum(xpsisquare) * x_unit**2 return end C_________EXPECTED Y-VALUE____________________ real function expected_y(psi,x_unit,length) implicit none integer i, j, length real x_unit complex, dimension(length,length):: psi,ypsisquare CMF$ LAYOUT psi(:news,:news) CMF$ LAYOUT ypsisquare(:news,:news) forall(i=1:length,j=1:length) & ypsisquare(i,j) = j * x_unit * abs(psi(i,j))**2 expected_y = sum(ypsisquare) * x_unit**2 return end C_______________________________________________________________________