Eigensystem calculates the eigen-energies and eigen-vectors of an Operator or matrix. The function comes in several flavors.

val, fun = Eigensystem(A)

• A : Square matrix (table of tables) of doubles

• val : Vector (Table of length #A) of doubles representing the eigen values
• fun : Matrix (Table of Tables of length #A by #A) representing the eigen vectors (fun[1] represents the first eigenvector)

A small example:

Eigensystem_Matrix.Quanty

A = {{ 1 , I , 0 },
     {-I , 2 , 1 },
     { 0 , 1 , 5 }}
 
val, fun = Eigensystem(A)
 
print(val)
print(fun)

Eigensystem_Matrix.out

{ 0.31866935639502 , 2.3579263675185 , 5.3234042760865 }
{ { (0 - 0.82050111444738 I) , 0.55903255238504 , -0.11941744665028 } , 
  { (0 - 0.56721932561261 I) , -0.77024207841542 , 0.29152937637548 } , 
  { (0 + 0.070994069063423 I) , 0.30693606176558 , 0.94907855109346 } }

psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)

calculates the lowest Npsi eigenstates for an operator Hamiltonian. The program uses iterative methods and needs to have a starting point to be defined. Note that the operator Hamiltonian does not know about the number of electrons there are in the system, only the number of orbitals is defined. If these orbitals are occupied or not is not given in the operator. A set of start restrictions is given by a table listing the number of Fermionic modes, Bosonic modes followed by lists including a string that tells the program which orbitals are included in the counting and then a minimum and maximum number of electrons in those orbitals. For example:

StartRestrictions = {Nf, Nb,  {"111111",2,2}}

Which would define two electrons in a p-shell. A start restriction of:

StartRestrictions = {Nf, Nb,  {"101010",2,2}}

would tell the program to start from two electrons with spin down and either 0,1,2, or 3 electron with spin up. A well defined calculation,but normally not a calculation one would like to do. If you want to start with two electrons with spin down and no electrons with spin up one needs:

StartRestrictions = {Nf, Nb,  {"101010",2,2}, {"010101",0,0}}

The final result must not fulfill the start restrictions, The program finds the lowest states in the basis defined by $(H+1)^{\infty} \psi_0$. Whereby $\psi_0$ is some random state that satisfies the start restrictions. If $S_z$ is a conserved quantum number of $H$ one would keep the number of spin up and down electrons fixed. If not other states with different values of $S_z$ will mix.

• Hamiltonian : Operator
• StartRestrictions : Table of restrictions (same format as restrictions that can be set for operators)
• Npsi : Number of lowest eigenstates to calculate. Calculation speed and memory uses depends on this number. States that are higher in energy then $\approx 10\times k_B \times T$ are not important for any physics and do not need to be calculated.

• psiList : A table of Wavefunctions of length Npsi

A small example:

Eigensystem_Startrestrictions.Quanty

dofile("../definitions.Quanty")
 
-- we define an Hamiltonian
Hamiltonian = OppU + 0.1 * Oppldots + 0.000001 * (2*OppSz+OppLz)
 
-- number of states calculated
Npsi=15
 
-- In order to make sure we have a filling of 2
-- electrons we need to define some restrictions
StartRestrictions = {Nf, Nb,  {"111111",2,2}}
 
-- We now can create the lowest Npsi eigenstates:
psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)
 
-- We define a list of some operators to look at their expectation values
oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz}
 
-- after we've created the eigen states we can look 
-- at a list of their expectation values
print(" <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>");
for i=1,#psiList do
  for j=1,#oppList do
    io.write(string.format("%7.4f ",psiList[i]*oppList[j]*psiList[i]))
  end
  io.write("\n")
end

Eigensystem_Startrestrictions.out

 <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>
-2.1033  1.9989  1.9989  0.0000 -0.0000  0.0000 
-2.0500  2.0000  2.0000  2.0000 -0.5000 -0.5000 
-2.0500  2.0000  2.0000  2.0000  0.0000 -0.0000 
-2.0500  2.0000  2.0000  2.0000  0.5000  0.5000 
-1.9521  1.9982  2.0036  6.0000 -0.9991 -1.0009 
-1.9521  1.9982  2.0036  6.0000 -0.4995 -0.5005 
-1.9521  1.9982  2.0036  6.0000  0.0000 -0.0000 
-1.9521  1.9982  2.0036  6.0000  0.4996  0.5004 
-1.9521  1.9982  2.0036  6.0000  0.9991  1.0009 
 0.4021  0.0018  5.9964  6.0000 -0.0009 -1.9991 
 0.4021  0.0018  5.9964  6.0000 -0.0005 -0.9995 
 0.4021  0.0018  5.9964  6.0000  0.0000 -0.0000 
 0.4021  0.0018  5.9964  6.0000  0.0005  0.9995 
 0.4021  0.0018  5.9964  6.0000  0.0009  1.9991 
 4.0033  0.0011  0.0011 -0.0000  0.0000 -0.0000 

Eigensystem(Hamiltonian, psiList)

Instead of starting from a random state that fulfills restrictions one can start from a pre-defined state. This can be useful if you have a large system and only want to do perturbation theory, configuration interaction for example, or if you want to calculate the eigen states of the Hamiltonian after making a small change to the Hamiltonian. In the example before we first calculate the eigen-state of the Hamiltonian without spin-orbit coupling and in a next step start from these eigen-states to calculate the Hamiltonian with spin-orbit coupling.

• Hamiltonian : Operator
• psiList : A table of Wavefunctions overwritten on output

• psiList : A table of Wavefunctions

A small example:

Eigensystem_Wavefunction.Quanty

dofile("../definitions.Quanty")
 
-- we define an Hamiltonian
Hamiltonian = OppU + 0.000001 * (2*OppSz+OppLz)
 
-- number of states calculated
Npsi=15
StartRestrictions = {Nf, Nb,  {"111111",2,2}}
psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)
 
oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz}
 
print(" <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>");
for i=1,#psiList do
  for j=1,#oppList do
    io.write(string.format("%7.4f ",psiList[i]*oppList[j]*psiList[i]))
  end
  io.write("\n")
end
 
Hamiltonian = Hamiltonian + 0.1 * Oppldots
 
oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz}
 
 
-- Recalculate eigen states with spin orbit
Eigensystem(Hamiltonian,psiList)
 
print(" <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>");
for i=1,#psiList do
  for j=1,#oppList do
    io.write(string.format("%7.4f ",psiList[i]*oppList[j]*psiList[i]))
  end
  io.write("\n")
end

Eigensystem_Wavefunction.out

 <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>
-2.0000  2.0000  2.0000  6.0000 -1.0000 -1.0000 
-2.0000  2.0000  2.0000  4.0000 -1.0000  0.0000 
-2.0000  2.0000  2.0000  4.0000  0.0000 -1.0000 
-2.0000  2.0000  2.0000  2.0000 -1.0000  1.0000 
-2.0000  2.0000  2.0000  4.0000  0.0000  0.0000 
-2.0000  2.0000  2.0000  2.0000  1.0000 -1.0000 
-2.0000  2.0000  2.0000  4.0000  0.0000  1.0000 
-2.0000  2.0000  2.0000  4.0000  1.0000  0.0000 
-2.0000  2.0000  2.0000  6.0000  1.0000  1.0000 
 0.4000  0.0000  6.0000  6.0000  0.0000 -2.0000 
 0.4000  0.0000  6.0000  6.0000  0.0000 -1.0000 
 0.4000  0.0000  6.0000  6.0000  0.0000  0.0000 
 0.4000  0.0000  6.0000  6.0000  0.0000  1.0000 
 0.4000  0.0000  6.0000  6.0000  0.0000  2.0000 
 4.0000 -0.0000 -0.0000 -0.0000  0.0000  0.0000 
 <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>
-2.1033  1.9989  1.9989  0.0000 -0.0000  0.0000 
-2.0500  2.0000  2.0000  2.0000 -0.5000 -0.5000 
-2.0500  2.0000  2.0000  2.0000  0.0000 -0.0000 
-2.0500  2.0000  2.0000  2.0000  0.5000  0.5000 
-1.9521  1.9982  2.0036  6.0000 -0.9991 -1.0009 
-1.9521  1.9982  2.0036  6.0000 -0.4995 -0.5005 
-1.9521  1.9982  2.0036  6.0000  0.0000 -0.0000 
-1.9521  1.9982  2.0036  6.0000  0.4996  0.5004 
-1.9521  1.9982  2.0036  6.0000  0.9991  1.0009 
 0.4021  0.0018  5.9964  6.0000 -0.0009 -1.9991 
 0.4021  0.0018  5.9964  6.0000 -0.0005 -0.9995 
 0.4021  0.0018  5.9964  6.0000  0.0000 -0.0000 
 0.4021  0.0018  5.9964  6.0000  0.0005  0.9995 
 0.4021  0.0018  5.9964  6.0000  0.0009  1.9991 
 4.0033  0.0011  0.0011 -0.0000  0.0000 -0.0000 

The last element of Eigensystem can be a table of options. Possible options are:

• “DenseBorder” Number of determinants in the bases where we switch from dense methods to sparse methods. (Standard value 1000)
• “NKrylovStart” Number of Krylov states in the Krylov basis in the first iteration. (Standard value 100)
• “NKrylovMax” Maximum number of Krylov states in the basis before a next iteration starts. (Standard value 500)
• “NKrylovStep” Increase in size of Krylov basis between each iteration. (Standard value 10)
• “NKrylovMin” Minimum number of states in the Krylov basis when running out of memory before the program stops storing the Krylov basis and goes to a mode where the states are recalculated (slower but much less memory). (Standard value 50)
• “NBitsKey” Number of bits in the key for the Hash lookup tables. (Standard value 16)
• “Epsilon” smallest coefficient before a determinant that is kept in the basis. (Standard value $1.49012*10^{-8}$)
• “Zero” Convergence limit (converged when variance of energy is smaller than epsilon). (Standard value $1.49012*10^{-6}$)
• “Restrictions” Possible set of restrictions for calculating the ground-state. (Standard value nil)
• “ExpandBasis” Boolean determining if the basis set should be kept fixed or allowed to expand. (Standard value true)