Yorick Quellcode Eigenfunktionen im Potentialgebirge

From Institute for Theoretical Physics II / University of Erlangen-Nuremberg

Jump to: navigation, search

Der folgende Quellcode kann in der frei verfügbaren Interpretersprache "Yorick" ausgeführt werden. Speichern Sie den Code z.B. als "eigenfunctions.i" ab und führen Sie ihn in Yorick aus, durch include, "eigenfunctions.i" .

Klicken Sie dann mit der Maus auf das Fenster, um das Potential zu verändern und automatisch die neuen Eigenfunktionen berechnen zu lassen!

Viel Spaß!

Florian Marquardt

// Yorick code for calculating the
// eigenenergies and eigenfunctions in an arbitrary
// 1D potential landscape. This works by
// setting up a tight-binding model (nearest neighbor hopping
// on a grid) and diagonalizing the resulting Hamiltonian matrix.
// The potential can be changed using the mouse.
// The program plots the local density of states (|Psi(x)|^2 at each eigenenergy)
// horizontal axis = x
// vertical axis = energy E
//
// by Florian Marquardt, 23.10.2011

// Parameters (numerical & physical)

N=200; // grid points in x-direction (= dimension of Hilbert space)
m=1.0; // mass (note hbar==1)

xmax=10.0; // extent of position interval (-xmax to xmax)

// plot the local density of states in the following energy interval:
Emin=-.1;
Emax=3.0;
// use this number of grid points for the energy direction (vertical axis):
NE=400;
// broaden the levels in energy by the following amount (to make them visible):
levelwidth=.01;


x=span(-xmax,xmax,N);

// preset functions

// particle in a box
func box {
  extern V, N, NE, levelwidth, Emin, Emax, xmax; 

  V=array(0.,N);
  Emin=-.1; Emax=3.; NE=400; levelwidth=.01;
}

// harmonic oscillator
func harmonic {
  extern V, N, NE, levelwidth, Emin, Emax, xmax; 

  omega=1.0;
  V=m*omega^2*x^2/2.;
  Emin=-.1; Emax=30.; NE=400; levelwidth=.1;
}

// linear rising potential
func linear {
  extern V, N, NE, levelwidth, Emin, Emax, xmax; 

  dVdx=1.0;
  V=dVdx*(x+xmax);
  Emin=0.; Emax=max(V); NE=400; levelwidth=.1;
}

// periodic
func periodic {
  extern V, N, NE, levelwidth, Emin, Emax, xmax; 

  Vamp=3.2;
  lambda=3.0;
  V=Vamp*(1+cos(x*2*pi/lambda));
  Emin=-.1; Emax=30.; NE=400; levelwidth=.1;
}

// just call those predefined functions here:
harmonic;

// an eigenvalue routine (using "singular value decomposition")
func SVeigen(a, &u)
/* DOCUMENT s = SVeigen(a, u)
*   return eigenvalues S and eignvectors U of symmetric matrix A.
*   That is,
*     A(,+) * U(+,) = S(-,) * U
*   The eigenvalues in S are in decreasing order of absolute value.
*   The eigenvectors are orthonormal:
*     U(+,) * U(+,) = unit(dimsof(A)(2))
*   If A is not symmetric, SVeigen returns garbage.
* SEE ALSO: SVdec
*/
{
  local vt;
  n = dimsof(a)(2);
  u = array(0., n, n);
  u(1:numberof(u):n+1) = x = 0.4964726826642538*avg(abs(a));
  /* this would work with x=0 except for case that A has
   * equal singular values corresponding to eigenvalues of opposite sign
   * adding essentially random x to all eigenvalues makes this very unlikely
   */
  s = SVdec(a+u, u, vt);          /* singular value = abs(eigenvalue) */
  i = abs(u+transpose(vt))(sum,) > abs(u)(sum,);    /* eigenvalue > 0 */
  s = (i+i-1.0)*s - x;
  /* sort unnecessary when x=0 */
  list = sort(-abs(s));
  s = s(list);
  u = u(,list);
  return s;
}



a=x(2)-x(1); // lattice spacing

// kinetic energy = hopping between adjacent sites
jump_matrix_element=1./(2*m*a^2);


// one run: generate Hamiltonian matrix and diagonalize it and plot result
func go {
  extern H, V, N, xmax, Emin, Emax, jump_matrix_element, Energies, DOS, x, X, E, U; 
  // now produce Hamiltonian:
  
  H=array(double,[2,N,N]);
  

  // we have set hbar=1
  for(j=1;j<N;j++) {
    H(j+1,j)=H(j,j+1)=-jump_matrix_element; // nearest neighbor hopping
    H(j,j)=2*jump_matrix_element; // diagonal elements
    H(j,j)+=V(j); // potential on diagonal entries
  }
  
  H(N,N)+=2*jump_matrix_element; 
  H(N,N)+=V(N);
  
  U=H; // matrix that will contain eigenvectors
  
// diagonalize the Hamiltonian, return eigenenergies and -vectors
  Energies=SVdec(H, U);
  
  
  // plot local density of states:
  fma;

  DOS=array(double,[2,N,NE]);
  // 2D array that will contain |Psi(x)|^2 at
  // different energies (=space-resolved or local density of states)
  
  E=span(Emin,Emax,NE)(-:1:N,); // energies (vertical direction)
  X=x(,-:1:NE); // x-coordinates (horizontal direction)
  
  // loop through all eigenvectors and add their |Psi|^2 to the plot
  for(j=1;j<=N;j++) {
    DOS+=exp(-((E-Energies(j))/levelwidth)^2)*U(,j)^2;
  }
  
  // now plot this
  fma; palette, "heat.gp";
  DOS=DOS^0.3;
  pli, DOS, -xmax, Emin, xmax, Emax;
  limits, -xmax, xmax, Emin, Emax;
  plg, V, x, color="white";
}



// calculate once:
go;

// now loop, always waiting for mouse clicks
// stop this infinite loop simply by interrupting the program (Ctrl-C)

editwidth=1.0;
stop=[];
do {
  result=mouse(-1, 0, ".");
  xpos=result(1);
  ypos=result(2);
  r=exp(-((x-xpos)/editwidth)^2);
  V=ypos*r+(1-r)*V;
  go;
 } while(stop==[]);