Differences

This shows you the differences between two versions of the page.

Link to this comparison view

documentation:tutorials:model_examples_from_physics:harmonic_oscillator [2016/10/07 21:15] – created Maurits W. Haverkortdocumentation:tutorials:model_examples_from_physics:harmonic_oscillator [2016/10/10 09:41] (current) – external edit 127.0.0.1
Line 1: Line 1:
 +{{indexmenu_n>2}}
 +====== Harmonic oscillator ======
  
 +###
 +This example solves the Harmonic oscillator numerically on a basis of plane waves.
 +###
 +
 +<code Quanty Harmonic_oscillator.Quanty>
 +-- solves the Harmonic oscillator H = -1/2 (d^2/dx^2) + 1/2 x^2
 +-- on a basis of complex plane waves
 +-- the plane wave basis assumes a periodicity, this length is:
 +a = 20
 +-- maximum k (ikmax * 2 * pi/a)
 +ikmax = 60
 +-- each plane wave is a basis "spin-orbital" k runs from -kmax to kmax, including 0, i.e. the number of basis "spin-orbitals" is:
 +NF = 2 * ikmax + 1
 +-- integration steps
 +dxint = 0.0001
 +
 +-- we first define a set of functions that are used to create the operators using integrals over the wave-functions
 +
 +-- the basis functions (plane waves) are:
 +function Psi(x, i)
 +  k = 2 * pi * i / a
 +  return (math.cos(k*x) + I * math.sin(k*x)) / math.sqrt(a)
 +end
 +
 +-- evaluate <psi_i | -1/2 (d^2/dx^2) | psi_j>
 +function IntegrateKineticEnergy(i,j)
 +  kj = 2 * pi * j / a
 +  sum = 0
 +  for x=-a/2, a/2, dxint do
 +    sum = sum - Conjugate(Psi(x,i)) * (kj*kj/2) * Psi(x,j) * dxint
 +  end
 +  return sum
 +end
 +
 +-- the previous integral has an analytical solution
 +function IntegrateKineticEnergyAna(i,j)
 +  if i==j then
 +    return 2*(j*pi/a)^2
 +  else
 +    return 0
 +  end
 +end
 +
 +-- evaluate <psi_i | 1/2 x^2 | psi_j>
 +function IntegratePotentialEnergy(i,j)
 +  kj = 2 * pi * j / a
 +  sum = 0
 +  for x=-a/2, a/2, dxint do
 +    sum = sum + Conjugate(Psi(x,i)) * (x*x/2) * Psi(x,j) * dxint
 +  end
 +  return sum
 +end
 +
 +-- the previous integral has an analytical solution
 +function IntegratePotentialEnergyAna(i,j)
 +  if i==j then
 +    return (a^2)/24
 +  else
 +    return a^2 * (2*(i-j) * pi * math.cos((i-j)*pi) + (-2+(i-j)^2*pi^2) * math.sin((i-j)*pi)) / (8 * (i-j)^3 * pi^3)
 +  end
 +end
 +
 +-- create the Hamiltonian
 +H = 0*NewOperator("CrAn",NF,0,0)
 +for i=-ikmax, ikmax do
 +  for j=-ikmax, ikmax do
 +    H = H + (IntegrateKineticEnergyAna(i,j) + IntegratePotentialEnergyAna(i,j)) * NewOperator("CrAn", NF, i+ikmax, j+ikmax)
 +  end
 +end
 +
 +-- create the operator to measure the kinetic energy
 +Hkin = 0*NewOperator("CrAn",NF,0,0)
 +for i=-ikmax, ikmax do
 +  for j=-ikmax, ikmax do
 +    Hkin = Hkin + IntegrateKineticEnergyAna(i,j) * NewOperator("CrAn", NF, i+ikmax, j+ikmax)
 +  end
 +end
 +
 +-- create the operator to measure the potential energy
 +Hpot = 0*NewOperator("CrAn",NF,0,0)
 +for i=-ikmax, ikmax do
 +  for j=-ikmax, ikmax do
 +    Hpot = Hpot + IntegratePotentialEnergyAna(i,j) * NewOperator("CrAn", NF, i+ikmax, j+ikmax)
 +  end
 +end
 +
 +-- The operators "do not know" how many electrons there are in the system.
 +-- We can restrict the solutions to only find states with one electron
 +-- For this we create a string with 1's and 0's that tells the code which basis orbitals to include for this specific restriction (useful if one has core states)
 +-- and in restrictions we set the minimum and maximum occupation of these basis orbitals (min=1 and max=1 in our case)
 +restrictionstring = ""
 +for i=1, NF do
 +  restrictionstring = restrictionstring .. "1"
 +end
 +StartRestrictions = {NF, 0, {restrictionstring,1,1}}
 +
 +
 +-- we now can create the lowest Npsi eigenstates:
 +Npsi=10
 +psiList = Eigensystem(H, StartRestrictions, Npsi)
 +-- note that for small systems Eigensystem creates a dense matrix and uses dense methods. 
 +-- The boarder between smal and big can be set, if you want to test the Lanczos algorithm as well
 +-- psiList = Eigensystem(H, StartRestrictions, Npsi,{{"DenseBoarder",5}})
 +
 +-- create output that lists the energy, kinetic energy and potential energy of the states found
 +oppList={H,Hkin,Hpot}
 +
 +print(" <E>            <Ekin>         <Epot>")
 +for key,psi in pairs(psiList) do
 +  expvalue = psi * oppList * psi
 +  for k,v in pairs(expvalue) do
 +    io.write(string.format("%14.7E ",v))
 +  end;
 +  io.write("\n")
 +end
 +</code>
 +
 +==== Result ====
 +The output is:
 +<file Quanty_Output Harmonic_oscillator.out>
 +Start of BlockGroundState. Converge 10 states to an energy with relative variance smaller than  1.490116119384766E-06
 +
 +Start of BlockOperatorPsiSerialRestricted
 +Outer loop   1, Number of Determinants:       121       121 last variance  2.926817131871151E+03
 +Start of BlockOperatorPsiSerialRestricted
 +Start of BlockGroundState. Converge 10 states to an energy with relative variance smaller than  1.490116119384766E-06
 +
 +Start of BlockOperatorPsiSerial
 + <E>            <Ekin>         <Epot>
 + 5.0000000E-01  2.5000000E-01  2.5000000E-01 
 + 1.5000000E+00  7.5000000E-01  7.5000000E-01 
 + 2.5000000E+00  1.2500000E+00  1.2500000E+00 
 + 3.5000000E+00  1.7500000E+00  1.7500000E+00 
 + 4.5000000E+00  2.2500000E+00  2.2500000E+00 
 + 5.5000000E+00  2.7500000E+00  2.7500000E+00 
 + 6.5000000E+00  3.2500000E+00  3.2500000E+00 
 + 7.5000000E+00  3.7500000E+00  3.7500000E+00 
 + 8.5000000E+00  4.2500000E+00  4.2500000E+00 
 + 9.5000000E+00  4.7500000E+00  4.7500000E+00 
 +</file>
 +
 +===== Table of contents =====
 +{{indexmenu>.#1|msort}}

You are here: Quanty - a quantum many body script language » Documentation » Tutorials » Model examples from physics » Harmonic oscillator

Print/export