ONE_PART_QHO

 !---------------------------------------------------------------

!Program name: ONE_PART_QHO

!Original code: Giannozzi,P., Ercolessi, F. and Gironcoli, S.

!Modified by Caio Lima Firme

!One particle quantum harmonic oscillator

!PARABOLIC POTENTIAL WELL

!Objective: to plot the wave function for each node

!---------------------------------------------------------------

!The dimensionless physical quantities:

!Potential energy: V(x)=0.5X^2   

! Energy: E=1/2+n

!---------------------------------------------------------------

! Use of Numerov method to find the wave function

! Function fn from Numerov equation is

! Fn = 1 + K2*h^2/12 , where K2=2m/((H/2pi)^2)*(E-V)

! H=Planck constant

!---------------------------------------------------------------

  integer :: grid, i, icl

  integer :: n, hn, ncross, j, iter

  double precision :: xmax, h, fh12, norm, djump, fac

  double precision :: Vupper, Vlower, e, k

  double precision :: x(0:300), y(0:300), V(0:300), fn(0:300)

  character (len=20) :: fileout

  character (len=1) :: A

!------------------------------------------------------

! Boundary conditions:

  y(0) = 0.0d0

  x(0) = 0.0d0

  xmax=10

  grid=300

! The number 300 refers to the grid everywhere in the code

  h =  xmax/grid

! Make explicit one term of the fn function from the Numerov method: fh12

  fh12=h*h/12.0d0

 !------------------------------------------------------

!Input data file

  print *, 'output file name is (include .dat at the end): '

  read *, fileout

!-------------------------------------------------------------

! Record in the computer memory each value of x(i) and V(i)

  do i = 0, grid

     x(i) = i * h

     V(i) = 0.5d0 * x(i)*x(i)

     write (*,*) i, 'X(i)=', x(i),'V(i)=', V(i)

  end do

!---------------------------------------------------------------

! Beginning of searching process

10 search: do

  !Input the number of nodes, n

20  print *, 'give the number of nodes, n (n=0,1,2,3,): '

  read *, n

  if (n < 0) then

  go to 20

  end if

!---------------------------------------------------------------

! set initial lower and upper bounds of the potential V(x)

  ! Set trial energy E

     Vupper=maxval (V(0:300))

     Vlower=minval (V(0:300))

     e = 0.5d0 * (Vlower + Vupper)

     iter = 100

Write (*,*) "Vupper=", Vupper, "e(initial)=", e

!------------------------------------------------------------------

 ! BEGINNING OF THE RECURSIVE PROCESS  (STEP=j)

  iterate: do j = 1, iter

! set up the fn-function used by the Numerov method

! fn=1+[2(V-E)]h^2/12 without m (mass) and H/2pi (Planck constant)

! determine the position of its last crossing, i.e. change of sign

     ! fn < 0 means classically allowed region

     ! fn > 0 means classically forbidden region

     fn(0)=fh12*(2.0d0*(V(0)-e))

     icl=-1

     do i=1,grid

        fn(i)=fh12*2.0d0*(V(i)-e)

        if ( fn(i) == 0.0d0) fn(i)=1.d-15

! store the index 'icl' where the last change of sign has been found

        if ( fn(i) /= sign(fn(i),fn(i-1)) ) icl=i

     end do

 ! fn as required by the Numerov method

     fn = 1.0d0 - fn

! determination of the wave-function in the first two points

 !--------------------------------------------------------------------

! Routine to determine whether n (number of nodes) is even or odd

  ! Check the program EVEN-ODD

     hn = n/2

     ! if n is even, the wave function is even

     ! if n is odd,  the wave function is odd

      if (2*hn == n) then

        ! even number of n: wavefunction is even

        y(0) = 1.0d0

        ! assume f(-1) = f(1)

        y(1) = 0.5d0*(12.0d0-10.0d0*fn(0))*y(0)/fn(1)

     else

        ! odd number of n: wavefunction is odd

        y(0) = 0.0d0

        y(1) = h

     end if

  !-------------------------------------------------

   ! outward integration and count number of crossings

     ncross=0

     do i =1,icl-1

  y(i+1)=((12.0d0-10.0d0*fn(i))*y(i)-fn(i-1)*y(i-1))/fn(i+1)

     if ( y(i) /= sign(y(i),y(i+1)) ) ncross=ncross+1

     end do

     fac = y(icl)

  !-------------------------------------------------

! Routine to determine whether n (number of nodes) is even or odd

  ! Check the program EVEN-ODD

     if (2*hn == n) then

        ! even number of n: no node in x=0

        ncross = 2*ncross

     else

        ! odd number of n: node in x=0

        ncross = 2*ncross+1

     end if

  !--------------------------------------------

     ! check number of crossings

     ! FIRST ITERATION PROCESS

 !THE NUMBER OF CROSSINGS HAS TO BE EQUAL TO THE NUMBER OF NODES

!Bisection method (see BISECTION3 code in chapter two)

     if ( iter > 1 ) then

        if (ncross /= n) then

   ! Incorrect number of crossings: adjust energy

if ( j == 1) &    print '("step         Energy              n     Discontinuity      Vupper     Vlower     icl")'        

write (*,5)  j, e, ncross, Vupper, Vlower, icl

5 format (i5,f25.15,i5,f35.15,f35.15,i40)

           if (ncross > n) then

      ! Too many crossings: current energy is too high,

      ! lower the upper bound

              Vupper = e

           else

       ! Too few crossings: current energy is too low,

          ! raise the lower bound

              Vlower = e

           end if

           ! New trial value:

           e = 0.5d0 * (Vupper+Vlower)

! go to beginning of do loop, don't perform inward integration

           cycle

         end if

      end if

!The equality ncross = n was reached

!END OF THE FIRST ITERATION PROCESS

!--------------------------------------------------------------------------

 ! Correct number of crossings: proceed to inward integration

  ! Determination of the wave-function in the last two points

     ! assuming y(grid+1) = 0

     y(grid) = h

   y(grid-1) = (12.0d0-10.0d0*fn(grid))*y(grid)/fn(grid-1)

     do i = grid-1,icl+1,-1

   y(i-1)=((12.0d0-10.0d0*fn(i))*y(i)-fn(i+1)*y(i+1))/fn(i-1)

     end do

! -------------------------------------------------------------------------

  ! rescale function to match at the classical turning point (icl)

     fac = fac/y(icl)

     Write (*,*) "fac=", fac

     y(icl:) = y(icl:)*fac

!-------------------------------------------------------------------------

     !Normalization of the wave function

     !DOT_PRODUCT : scalar product of vectors

     !See the code: NORMALIZATION

     norm = dot_product (y, y)

     y = y / sqrt(norm)

!----------------------------------------------------------------------

! SECOND ITERATION PROCESS

! Bisection method (See BISECTION3 code in chapter one)

     if ( iter > 1 ) then

      ! calculate the discontinuity in the first derivative

        !  y'(i;RIGHT) - y'(i;LEFT)

 djump = ( y(icl+1) + y(icl-1) - (14.0d0-12.0d0*fn(icl))*y(icl) ) / h

     write (*,30) j, e, n, djump, Vupper, Vlower,icl

     30 format (i5,f25.15,i5,f14.8,f22.15,f34.15,i40)

          if (djump*y(icl) > 0.0d0) then

    ! Energy is too high --> choose lower energy range

           Vupper = e

        else

    ! Energy is too low  --> choose upper energy range

           Vlower = e

        endif

        e = 0.5d0 * (Vupper+Vlower)

! -------------------- convergence test -----------------------------

        if ( Vupper-Vlower < 1.d-10) exit iterate

     endif

  !END OF THE SECOND ITERATION PROCESS

  end do iterate

  !END OF THE RECURSIVE PROCESS

 !-------------Convergence achieved-----------------------

  print *, 'Do you want to do another calc (Y/N)?'

           read *, A

           Select case (A)

                 Case ('Y','y')

                   go to 10

                 case default

                   continue

           end select

!------------------------------------------------------------------

!Printing results  x(i) and y(i) in the 'fileout'

! x<0 region:

open(7,file=fileout,status='replace')

  do i=grid,1,-1

     write (7,'(f7.3,3e16.8,f12.6)')  -x(i), ((-1)**n)*y(i)

  enddo

  ! x>0 region:

  do i=0,grid

     write (7,'(f7.3,3e16.8,f12.6)') x(i), y(i)

  enddo

  write (7,'(/)')

  close(7)

  stop

!--------------------------------------------------------------------------

end do search

 

stop

end

Comentários

Postagens mais visitadas deste blog

RECURSIVE

HYDROGEN-RADIAL

RUNGE-KUTTA