! Time independent quantum mechanics - Lennard-Jones potential in one dimension
! solve using variational-monte carlo method
! Program to accompany "Computational Physics" by N. Giordano
! Copyright Prentice Hall 1997
program monte
   option nolet
   dim psi(0 to 2000)   ! wave function
   call initialize(psi,dx,d_psi,x_left,x_right,sigma,epsilon,energy)
   n_total = 0
   n_moves = 0
!   call save(psi,x_left,x_right,dx,energy,n_total)
   do
      call calculate(psi,dx,d_psi,x_left,x_right,sigma,epsilon,energy,n_accepted,n_total,n_moves)
      print energy,d_psi,n_accepted;n_total;n_moves
      call display(psi,dx,x_left,x_right)
      if key input then
         get key z
         c$ = chr$(z)
                           ! save the results if desired
         if c$ = "s" then call save(psi,x_left,x_right,dx,energy,n_total)
         if c$ = "c" then  ! clear screen if too cluttered
            clear
            call axes(sigma,epsilon)   ! plot potential for comparison
         end if
      end if
   loop
end
! initialize variables
! psi(i) = wave function at grid site i
! dx = grid step   d_psi = max size of Monte-Carlo changes in psi
! x_left,x_right = boundaries of region (psi = 0 outside this)
! sigma,epsilon = Lennard-Jones parameters
! energy = energy of initial trial function
sub initialize(psi(),dx,d_psi,x_left,x_right,sigma,epsilon,energy)
   dx = 0.2
   sigma = 1
   epsilon = 10 
   x_left = 0.5 * sigma
   x_right = 5 * sigma
   for i = x_left/dx to x_right/dx     ! initialize wave function
      psi(i) = 0
   next i
   delta = 3 * sigma / dx
   i_start = 1.1 * sigma / dx 
   for i = i_start to i_start + delta
      psi(i) = 1 / sqr(delta)
   next i
   d_psi = 0.1 / sqr(delta)
   call normalize(psi,dx,x_left,x_right)   ! normalize trial function
   call calc_energy(psi,dx,x_left,x_right,sigma,epsilon,energy)
   set window 0,6,-2,2              ! calculate energy of trial function
   call axes(sigma,epsilon)
   call display(psi,dx,x_left,x_right)  ! display wave function
end sub 
! plot axes and potential
sub axes(sigma,epsilon)
   declare def potential
   set color "black"
   clear
   max = 5
   plot 0,0;max,0
   for x = 0.7 to max step 0.1
      plot x,potential(x,sigma,epsilon) / 6;
   next x 
   plot x,potential(x,sigma,epsilon) / 6
end sub
! do the work here
! n_total = total number of Monte Carlo moves attempted
! n_moves = total number of Monte Carlo moves accepted
sub calculate(psi(),dx,d_psi,x_left,x_right,sigma,epsilon,energy,n_a,n_total,n_moves)
   declare def potential    ! declare my own function for this
   max = x_right/dx
   min = x_left/dx
   n_a = 0
   n_attempts = 1000        ! make this many attempted moves
   for i = 1 to n_attempts
      n = min + int(rnd*(max + 1 - min))  ! choose site to adjust randomly
      psi_old = psi(n)                    ! with uniform probability
      psi(n) = psi(n) + 2*(rnd-0.5)*d_psi ! adjust trial function
      call calc_energy(psi,dx,x_left,x_right,sigma,epsilon,new_energy)
      if new_energy > energy then     ! compare new energy to old energy
         psi(n) = psi_old             ! keep new trial function only if it
      else                            ! lowers the energy
         energy = new_energy
         call normalize(psi,dx,x_left,x_right)
         n_a = n_a + 1
      end if
   next i
   n_total = n_total + n_attempts
   n_moves = n_moves + n_a
end sub
! display wave function
sub display(psi(),dx,x_left,x_right)
   set color "red"
   for n = x_left/dx - 2 to x_right/dx + 2
      plot n*dx,psi(n)
   next n
   plot n*dx,psi(n)
end sub
! Lennard-Jones potential
def potential(x,sigma,epsilon)
   r6 = (x/sigma)^(-6)
   potential = 4 * epsilon * (r6^2 - r6)
end def  
! save wave function to a file
sub save(psi(),x_left,x_right,dx,e,n_tot)
   open #1: name "lj_monte." & str$(e) & "." & str$(n_tot), create new, organization text
   m_left = x_left / dx
   m_right = x_right / dx
   for i = m_left to m_right
      print #1: i*dx,psi(i)
   next i
   close #1
end sub
! calculate energy of trial function
sub calc_energy(psi(),dx,x_left,x_right,sigma,epsilon,energy)
   declare def potential
   energy = 0
   sum = 0
   for i = x_left/dx-1 to x_right/dx+1
      energy = energy + dx * potential(i*dx,sigma,epsilon)*psi(i)^2 - dx * (0.5/(dx^2))*psi(i)*(psi(i+1) + psi(i-1) - 2*psi(i))
      sum = sum + psi(i) * psi(i) * dx
   next i
   energy = energy / sum
end sub
! normalize a trial function
sub normalize(psi(),dx,x_left,x_right)
   sum = 0
   for i = x_left/dx-1 to x_right/dx+1
      sum = sum + dx * psi(i) * psi(i)
   next i
   for i = x_left/dx-1 to x_right/dx+1
      psi(i) = psi(i) / sqr(sum)
   next i
end sub