!=========================================================== ! ! File: sch.f90 ! ! Integrate 1d Schrodinger equation from xmin to xmax. ! Determine energy eigenvalue and eigenfunction by matching ! evolving solutions from xmin and from xmax at a point xm. ! Matching done by equating values of functions and their ! derivatives at xm. The point xm chosen at the left most ! turning point of the potential at any given value of the ! energy. The potential and boundary conditions chosen in ! different file. ! ---------------------------------------------------------- ! Input: energy: Trial value of energy ! de: energy step, if matching fails de -> e+de, if ! logderivative changes sign de -> -de/2 ! xmin, xmax, Nx ! ---------------------------------------------------------- ! Output: Final value of energy, number of nodes of ! wavefunction in stdout ! Final eigenfunction in file psi.dat ! All trial functions and energies in file all.dat ! ---------------------------------------------------------- ! Compile: ! gfortran schInfSq.f90 sch_util.f90 sch.f90 -o s !=========================================================== program schroedinger_equation_1D use, intrinsic :: iso_fortran_env use :: schV use :: schutil implicit none integer , parameter :: dp = real64 ! ---------------------------------------------------------- real(dp), allocatable :: psi(:),psip(:) integer :: Nx,NxL,NxR real(dp) :: dx real(dp) :: xmin,xmax,xm !left/right/matching points real(dp) :: psixmin,psipxmin,psixmax,psipxmax real(dp) :: psileft ,psiright ,psistep,psinorm real(dp) :: psipleft,psipright,psipstep real(dp) :: de,epsilon real(dp) :: matchlogd,matchold,psiold,norm,x integer :: iter,i,imatch,nodes integer :: uout !file unit for output !---------- Input: print *,'# Enter energy,de,xmin,xmax,Nx' read *, energy,de,xmin,xmax,Nx !--- need even intervals for normalization integration if( xmin >= xmax) stop 'Error: xmin >= xmax' if( mod(Nx,2) == 0 ) Nx=Nx+1 allocate(psi(Nx),psip(Nx)) dx = (xmax - xmin)/(Nx-1) epsilon = 1.0e-6_dp call boundary(xmin,xmax,psixmin,psipxmin,psixmax,psipxmax) print *,'# #######################################' print *,'# Estart= ',energy, ' de= ',de print *,'# Nx= ',Nx ,' eps= ',epsilon print *,'# xmin= ',xmin,' xmax= ',xmax, ' dx= ',dx print *,'# psi(xmin)= ',psixmin,' psip(xmin)= ',psipxmin print *,'# psi(xmax)= ',psixmax,' psip(xmax)= ',psipxmax print *,'# #######################################' !----- Calculate: open(newunit=uout,file='all.dat') matchold = 0.0_dp do iter = 1, 10000 !----- Determine matching point at turning point from the left: imatch = -1 do i = 1, Nx x = xmin + (i-1)*dx if( imatch < 0 .and. (energy-V(x)) > 0.0_dp) imatch = i enddo if ( imatch <= 100 .or. imatch >= Nx-100) imatch = Nx/5 xm = xmin + (imatch-1)*dx NxL = imatch NxR = Nx-imatch+1 !----- Evolve wavefunction from the left: psi (1) = psixmin psip (1) = psipxmin psistep = psixmin psipstep = psipxmin do i=2,NxL x = xmin + (i-2)*dx !this is x before the step call RKSTEP(x,psistep,psipstep, dx) psi (i) = psistep psip(i) = psipstep enddo ! use this to normalize eigenfunction to match at xm psinorm = psistep psipleft = psipstep !----- Evolve wavefunction from the right: psi (Nx) = psixmax psip(Nx) = psipxmax psistep = psixmax psipstep = psipxmax do i=2,NxR x = xmax - (i-2)*dx call RKSTEP(x,psistep,psipstep,-dx) psi (Nx-i+1) = psistep psip(Nx-i+1) = psipstep enddo psinorm = psistep/psinorm psipright = psipstep !----- Renormalize psil so that psil(xm)=psir(xm) psi (1:NxL-1) = psinorm * psi (1:NxL-1) psip(1:NxL-1) = psinorm * psip(1:NxL-1) psipleft = psinorm * psipleft !----- print current solution: do i = 1, Nx x = xmin + (i-1)*dx write(uout,*)iter,energy,x,psi(i),psip(i) enddo !----- matching using derivatives: !Careful: this can fail if psi'(xm) = 0 !! (use also |de|<1e-6 !criterion) matchlogd = & (psipright-psipleft)/(ABS(psipright)+ABS(psipleft)) print *,'# iter,energy,de,xm,logd: ',& iter,energy,de,xm,matchlogd !----- Exit condition: if(ABS(matchlogd) <= epsilon.or.ABS(de/energy) < 1.0e-12_dp)& exit if( matchlogd * matchold < 0.0_dp) de = -0.5_dp*de energy = energy + de matchold = matchlogd enddo ! do iter=1,10000 close(uout) !--------------------------------------------------- !----- Solution has been found and now it is stored: norm = integrate(psi,dx,Nx) norm = 1.0D0/sqrt(norm) psi(:) = norm*psi(:) !----- Cound number of zeroes, add one and get energy level: nodes = 1 psiold = psi(1) do i=2,Nx-1 if( ABS( psi(i)) > epsilon)then !should be 0 within epsilon if( psiold*psi(i) < 0.0_dp )nodes = nodes+1 psiold = psi(i) endif enddo !i=2,Nx-1 !------- Print final solution: open(newunit=uout,file='psi.dat') print *,'Final result: E= ',energy,' n= ',nodes,& ' norm = ' ,norm if( ABS(matchlogd) .gt. epsilon) print *& ,'Final result: SOS: logd>epsilon. logd= ',matchlogd write (uout,*)'# E= ' ,energy,' n= ',nodes,& ' norm = ' ,norm do i = 1, Nx x = xmin + (i-1)*dx write(uout,*) x,psi(i) enddo close (uout) end program schroedinger_equation_1D ! --------------------------------------------------------------------- ! Copyright by Konstantinos N. Anagnostopoulos (2004-2021) ! Physics Dept., National Technical University, ! konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant ! ! This program is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation, version 3 of the License. ! ! This program is distributed in the hope that it will be useful, but ! WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! General Public License for more details. ! ! You should have received a copy of the GNU General Public Liense along ! with this program. If not, see . ! -----------------------------------------------------------------------