
BeginPackage["SACI`", {"DiscreteMath`Combinatorica`", "LinearAlgebra`Orthogonalization`",
     "LinearAlgebra`MatrixManipulation`"}]
(*
SACI`
*)
NElectronCreate::usage  := "NElectronCreate[N eee] "

ShowSpinEigenfunctionCompatibility::usage = "Prints a table showing where a non zero element in row i column j means that the jth spin eigenfunction is compatiable with i doubly occupied orbitals."
\
(*
Prints a table showing where a non zero element in row i column j means that the jth spin eigenfunction is compatiable with i doubly occupied orbitals.
\
*)
energyconversionfactor::usage := "energyconversionfactor[relmass,relperm] gives the conversion factor between energy in material units and eV"
\!\(OneElectron::usage\  := \ \*"\"\<OneElectron[EE_] produces all pairs of one electron quantum numbers that have a zero magnetic field energy equal to EE*\[HBar] \!\(\[Omega]\_0\) \>\""\)
SingleOrbitalEnergyMagnetic::usage := "SingleOrbitalEnergyMagnetic[u_List, firengap_, magfield_, relmass_, relperm_]"

LoadIntegralandHamiltonianData::usage   := "Loads the intermediate steps in the calculation of the interaction intergral as well as the Interaction Hamiltonian matrix elements"

SaveIntegralandHamiltonianData::usage   := "Saves the intermediate steps in the calculation of the interaction intergral as well as the Interaction Hamiltonian matrix elements"

LevelVsField::usage =  "LevelVsField[Hint_, basis_, firengap_, relmass_,relperm_, Lvalue_ ,magmin_, magmax_, step_ ,levlist_] produces a list ofthe energy levels for the range of energies specified by magmin, magmax and step in Tesla."
\
(*
LevelVsField[Hint_, basis_, firengap_, relmass_,relperm_, Lvalue_ ,magmin_, magmax_, step_ ,levlist_] produces a list ofthe energy levels for the range of energies specified by magmin, magmax and step in Tesla.
\
*)
lengthscale::usage = "lengthscale[relmass_, relperm_] gives the effective Bohr radius from the relative effective mass and relativepermittivity."
\
(*
lengthscale[relmass_, relperm_] gives the effective Bohr radius from the relative effective mass and relativepermittivity.
\
*)
ElectronDensityFunction::usage = "ElectronDensityFunction[vec_, basisgen_,r_,\[Alpha]_, firstenergygap_, magneticfield_, relativeeffectivemass_, relativepermittivity_] calculates the electron probablity density as from the normalised energy eigenvector vec."
\
(*
ElectronDensityFunction[vec_, basisgen_,r_,\[Alpha]_, firstenergygap_, magneticfield_, relativeeffectivemass_, relativepermittivity_] calculates the electron probablity density as from the normalised energy eigenvector vec.
\
*)
Convergence::usage = 
  "Convergence[Hint_, basis_, firengap_, magfield_, relmass_,relperm_, Lvalue_, sizelist_List, levlist_] uses the Interaction Hamiltonian to produce a plot of the convergence for various basis sizes given by sizelist and for the level given by levlist."
\
(*
Convergence[Hint_, basis_, firengap_, magfield_, relmass_,relperm_, Lvalue_, sizelist_List, levlist_] uses the Interaction Hamiltonian to produce a plot of the convergence for various basis sizes given by sizelist and for the level given by levlist.
\
*)
EnergyLevelPlot::usage = "EnergyLevelPlot[Hfin_SparseArray, numberoflevels_,relmass_, relperm_,  opts___] plot the energy of the number of levels given"
\
(*
EnergyLevelPlot[Hfin_SparseArray, numberoflevels_,relmass_, relperm_,  opts___] plot the energy of the number of levels given
\
*)
TotalHamiltonian::usage = "TotalHamiltonian[Hint_SparseArray, basisgen_List, firengap_, magfield_, relmass_, relperm_] calculates the total Hamiltonian given the interaction Hamiltonian, the basis for the interaction Hamiltonian and the semiconductor and harmonic well parameters."
\
(*
TotalHamiltonian[Hint_SparseArray, basisgen_List, firengap_, magfield_, relmass_, relperm_] calculates the total Hamiltonian given the interaction Hamiltonian, the basis for the interaction Hamiltonian and the semiconductor and harmonic well parameters.
\
*)
InteractionHamiltonian::usage = "InteractionHamiltonian[basisgen_List] calculates the interaction Hamiltonian for a particular basis."
\
(*
InteractionHamiltonian[basisgen_List] calculates the interaction Hamiltonian for a particular basis.
\
*)
InteractionHamiltonianMatrixElement::usage = "InteractionHamiltonianMatrixElement[repa_AngleBracket, repb_AngleBracket] calculates the Hamiltonian matrix element between two basis elements."
\
(*
InteractionHamiltonianMatrixElement[repa_AngleBracket, repb_AngleBracket] calculates the Hamiltonian matrix element between two basis elements.
\
*)
InteractionIntegral::usage = 
  "InteractionIntegral[n,m,p,q,a,b,c,d] calculates the pairwise interaction integral for electrons obeying Fock-Darwin solutions interacting through a Coulomb potential for the quantum numbers n, m, p, q, a, b, c and d."
\
(*
InteractionIntegral[n,m,p,q,a,b,c,d] calculates the pairwise interaction integral for electrons obeying Fock-Darwin solutions interacting through a Coulomb potential for the quantum numbers n, m, p, q, a, b, c and d.
\
*)
ElectronList::usage = "The variable \!\(\*
StyleBox[\"ElectronList\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\" \",\nFontVariations->{\"CompatibilityType\"->0}]\)is a list where \!\(TraditionalForm\`1\) stands for a spin up electron and 0 stands for a spin down electron.  This implicitly defines the number of electrons by the length of the list and the total projection of spin angular momentum quantum number by the number of 1's and 0's. e.g. {1, 1, 0}";\

(*
The variable ElectronList is a list where "1" stands for a spin up electron \
and "0" stands for a spin down electron.  This implicitly defines the number \
of electrons by the length of the list and the total projection of spin \
angular momentum quantum number by the number of 1's and 0's.
*)
S::usage = "The square of angular momentum quantum number is given by \!\(\*
StyleBox[\"S\",\nFontWeight->\"Bold\"]\).  To get meaningful results, \!\(\*
StyleBox[\"S\",\nFontWeight->\"Bold\"]\) must be compatible with the number and spin of the electrons defined by  \!\(\*
StyleBox[\"ElectronList\",\nFontWeight->\"Bold\"]\).";
(*
The square of angular momentum quantum number is given by S.  To get \
meaningful results, S must be compatible with the number and spin of the \
electrons defined by  ElectronList.
*)
CalculateSpinBasis::usage = "CalculateSpinBasis[] calculates the spin eigenfunctions that form the spin eigenspace given the number of electrons and M quantum number define by ElectronList as well as the total spin Quantum Number S.";
  
SingleElectronProbablityDensity::usage = "SingleElectronProbablityDensity[n,m, firstenergygap_, magneticfield_, rmax_] plots the probability density of the Fock-Darwin solution with quantum numbers n and m with the harmonic well constant implicitly defined by firstenergygap in meV and a magnetic field in Tesla. rmax gives the radial range of the plot in multiples of the effective Bohr radius. relativeeffectivemass and relativepermittivity must be also be pre-defined.";
  
ShowSpinBasis::usage = "\!\(\*
StyleBox[\"ShowSpinBasis\",\nFontWeight->\"Bold\"]\) can display the spin basis after it has been calculated by \!\(\*
StyleBox[\"CalculateSpinBasis\",\nFontWeight->\"Bold\"]\)\!\(\*
StyleBox[\".\",\nFontWeight->\"Bold\"]\)";
RepresentationMatrix::usage  = "RepresentationMatrix[perm_List] displays the representation matrix associated associated with the permutation 'perm'.  For N electrons, the permutation is given by a reording of the list {1,2,3,..N}, e.g. {2,1,3,4,5} for the transposition of electrons 1 and 2 in a 5 electron system";
  
BasisCreate::usage = "Creates a basis for Matrix diagonalisation";
BasisCreateOneL::usage = \
"Creates a basis for Matrix diagonalisation with only one sum of orbital angular momentum quantum number";\

\[Alpha]::usage  = "Denotes spin up electron."
(*
Denotes spin up electron.
*)
\[Beta]::usage = "Denotes spin down electron."
(*
Denotes spin down electron.
*)
NElectronDouble::usage = ".";
a::usage = "";
b::usage = "";
(*
3.2 Basis Creation Functions
*)
\!\(\(\(OneElectron[EE_]\)\(\ \)\(:=\)\(\ \)\(If[OddQ[EE], \ Delete[Flatten[Table[{n, \ i \((EE - 1 - 2  n)\)}, \ {n, 0, \(EE - 1\)\/2\ }, {i, \(-1\), 1, 2}], 1], EE + 1], \ Flatten[Table[{n, i \((EE - 1 - 2  n)\)}, {n, 0, \(EE - 2\)\/2\ }, {i, \(-1\), 1, 2}], 1]]\)\(\ \)\)\)
(*
OneElectron[EE_] produces a list of all pairs of n,m quantum numbers that \
give a one electron energy of EE.
*)
OneElectronBasis[EE_] := Flatten[Table[OneElectron[ee], {ee, 1, EE}], 1]
(*
OneElectronBasis[EE_] produces a list of all pairs of n,m quantum numbers \
that give a one electron energy less than or equal to EE.
*)
OneElectronEnergyList[E_] := Flatten[Table[e, {e, 1, E}, {n, 1, e}]]
(*
OneElectronEnergyList[EE_] produces a list of energies that corresponds to \
the list of quantum number pairs given by OneElectronBasis[EE_]  for the same \
value of EE.
*)
OneElectronBasisSize[EE_, NN_, D_] := 
  EE - Plus @@ (Sort[
            OneElectronEnergyList[NN] ~Join~ 
              OneElectronEnergyList[NN][[Range[1, D]]]])[[Range[1, NN - 1]]] 
NElectronSingle[NN_, EE_, fff_] := 
  Module[{String1, String2, String3, String4, i, rule1}, 
    String1 = "";
    For[i = 1, i <= NN, 
      String1 = 
        StringJoin[String1, 
          StringReplace[
            "___, {ri_,rj_}, ", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String1 = StringJoin["{", String1, "___}"];
    String2 = "";
    For[i = 1, i <= NN - 1, 
      String2 = 
        StringJoin[String2, 
          StringReplace[
            "{ri,rj},", {"i" -> ToString[2i - 1], "j" -> ToString[2i] }]]; 
      i++];
    String2 = 
      StringJoin["\[LeftAngleBracket]{", String2, 
        StringReplace[
          "{ri,rj}", {"i" -> ToString[2 NN - 1], "j" -> ToString[2 NN] }], 
        "},k\[RightAngleBracket]"];
    String3 = "";
    For[i = 1, i <= fff, 
      String3 = 
        StringJoin[String3, 
          StringReplace[If[i == fff, String2, StringJoin[String2, ","]], 
            "k" -> ToString[i]]]; i++];
    String3 = StringJoin["{", String3, "}"];
    String4 = "";
    For[i = 1, i <= NN - 1, 
      String4 = 
        StringJoin[String4, 
          StringReplace[
            "2*ri + Abs[rj] + 1 + ", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String4 = 
      StringJoin["(", String4,  
        StringReplace[
          "2*ri + Abs[rj] + 1", {"i" -> ToString[2 NN - 1], 
            "j" -> ToString[2 NN] }], 
        StringReplace[" == ee)", "ee" -> ToString[EE]]];
    rule1 = ToExpression[ StringJoin[String1, "/;", String4, "->", String3]];
    Flatten[
      ReplaceList[OneElectronBasis[OneElectronBasisSize[EE, NN, 0]], rule1], 
      1]]
OneElectronPart[NN_, EE_] := 
  Module[{String1, String2, String3, String4, i, rule1},  
    String1 = "";
    For[i = 1, i <= NN, 
      String1 = 
        StringJoin[String1, 
          StringReplace[
            "___, {ri_,rj_}, ", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String1 = StringJoin["{", String1, "___}"];
    String2 = "";
    For[i = 1, i <= NN - 1, 
      String2 = 
        StringJoin[String2, 
          StringReplace[
            "{ri,rj},", {"i" -> ToString[2i - 1], "j" -> ToString[2i] }]]; 
      i++];
    String2 = 
      StringJoin["\[LeftAngleBracket]", String2,  
        StringReplace[
          "{ri,rj}", {"i" -> ToString[2 NN - 1], "j" -> ToString[2 NN] }], 
        "\[RightAngleBracket]"];
    String4 = "";
    For[i = 1, i <= NN - 1, 
      String4 = 
        StringJoin[String4, 
          StringReplace[
            "2*ri + Abs[rj] + 1 + ", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String4 = 
      StringJoin["(", String4,  
        StringReplace[
          "2*ri + Abs[rj] + 1", {"i" -> ToString[2 NN - 1], 
            "j" -> ToString[2 NN] }], 
        StringReplace[" == ee)", "ee" -> ToString[EE]]];
    rule1 = ToExpression[ StringJoin[String1, "/;", String4, "->", String2]];
    ReplaceList[OneElectronBasis[OneElectronBasisSize[EE, NN, 0]], rule1]]
DoubleElectronPart[DD_, EE_] := 
  Module[{String1, String2, String3, String4, i, rule1},  
    String1 = "";
    For[i = 1, i <= DD, 
      String1 = 
        StringJoin[String1, 
          StringReplace[
            "___, {ri_,rj_}, ", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String1 = StringJoin["{", String1, "___}"];
    String2 = "";
    For[i = 1, i <= DD - 1, 
      String2 = 
        StringJoin[String2, 
          StringReplace[
            "{ri,rj},{ri,rj},", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String2 = 
      StringJoin["\[LeftAngleBracket]", String2,  
        StringReplace[
          "{ri,rj},{ri,rj}", {"i" -> ToString[2 DD - 1], 
            "j" -> ToString[2 DD] }], "\[RightAngleBracket]"];
    String4 = "";
    For[i = 1, i <= DD - 1, 
      String4 = 
        StringJoin[String4, 
          StringReplace[
            "2(2*ri + Abs[rj] + 1) + ", {"i" -> ToString[2i - 1], 
              "j" -> ToString[2i] }]]; i++];
    String4 = 
      StringJoin["(", String4,  
        StringReplace[
          "2(2*ri + Abs[rj] + 1)", {"i" -> ToString[2 DD - 1], 
            "j" -> ToString[2 DD] }], 
        StringReplace[" == ee)", "ee" -> ToString[EE]]];
    rule1 = ToExpression[ StringJoin[String1, "/;", String4, "->", String2]];
    ReplaceList[OneElectronBasis[OneElectronBasisSize[EE, 2DD, DD]], rule1]]
Trial[dd_, ee_, EE_] := 
  DeleteCases[
    Flatten[Outer[Join, DoubleElectronPart[dd, ee], 
        OneElectronPart[NN - 2dd, EE - ee], 1], 
      1], \[LeftAngleBracket]___, {a_Integer, b_Integer}, ___, {a_Integer, 
        b_Integer}, ___, {a_Integer, b_Integer}, ___\[RightAngleBracket], 1]
NElectronDouble[NN_, EE_ , dd_, ff_] := Module[{}, ZZ = {};
    If[NN - 2dd != 0,
      
      For[ee = 2, ee <= EE, 
          ZZ = Join[ZZ, 
              DeleteCases[
                  Flatten[
                    Outer[Join, DoubleElectronPart[dd, ee], 
                      OneElectronPart[NN - 2dd, EE - ee], 1], 
                    1], \[LeftAngleBracket]___, {a_, b_}, ___, {a_, 
                      b_}, ___, {a_, b_}, ___\[RightAngleBracket], 
                  1] /. \[LeftAngleBracket]x___\[RightAngleBracket] -> \
\[LeftAngleBracket]{x}, ff\[RightAngleBracket] ]; ee = ee + 2 ];,
      ZZ = 
        Join[ZZ, 
          DoubleElectronPart[dd, 
              EE] /. \[LeftAngleBracket]x___\[RightAngleBracket] -> \
\[LeftAngleBracket]{x}, ff\[RightAngleBracket]]];  ZZ]
NElectronCreate[numelectrons_, energy_] := 
  Module[{basis, icount, jjcount}, 
    basis = NElectronSingle[numelectrons, energy, Length[Vtrans]]; 
    For[icount = 1, icount <= numpairs, 
      For[jjcount = 1, jjcount <= Length[Vtrans],
        If[Norsym[[icount, jjcount]] !=  0, 
          basis = Join[basis, 
              NElectronDouble[numelectrons, energy, icount, jjcount ] ] ] ; 
        jjcount++]; icount++]; basis]
NElectronBasisCreate[Numelec_, Ene_] := 
  Module[{Nbasis, eee}, NBasis = {}; 
    For[eee = Numelec, eee <= Ene, 
      NBasis = Join[NBasis, NElectronCreate[Numelec, eee]]; eee++]; NBasis]
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(Nbasis\\)\\\" is similar to existing symbol \\\"\\!\\(basis\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(NBasis\\)\\\" is similar to existing symbol \\\"\\!\\(Nbasis\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
\!\(NElectronBasisCreateOneL[Numelec_, \ Ene_, \ Lvalue_]\  := Select[\ NElectronBasisCreate[Numelec, Ene], \((\[Sum]\+\(i = 1\)\%NN\((#1[\([1, i, 2]\)])\)\  == \ \ Lvalue)\) &]\)
BasisCreate[Ene_] := NElectronBasisCreate[NN, Ene]
BasisCreateOneL[Ene_, Lvalue_] := NElectronBasisCreateOneL[NN, Ene, Lvalue]
(*
In the non-interacting electron case, the multiple electron energy \
eigenstates are simply the antisymmetrized products of the one-electron \
eigenstates with spin-eigenfunctions. The energy eigenvalue is then the sum \
of the one electron energies. In the interacting case, the energy eigenstates \
can be approximated as a finite linear combination of the non-interacting \
eigenstates with the approximation improving as the number of non-interacting \
eigenstates included. 
The above functions are used to calculate a basis consisting of all possible \
spatial and spin combinations with a total non-interacting energy less than \
or equal to a certain value. This proves to be an effective basis with \
convergence improving as the cut off energy is increased.
*)
(*
Begin
*)
Begin["`Private`"]
(*
SACI`Private`
*)
(*
Misc, Data
*)
(*
Messages
*)
(*
This turns off some error messages.
*)
Off[General::"spell"]; Off [General::"spell1"];
Off[Unset::"norep"];
Off[Interpolation::"inhr"];
(*
Physical Constants
*)
(*
The following define the physical constants used in the program.
*)
\!\(\(electroncharge\  = \ 1.602176462*10\^\(-19\);\)\)
\!\(\(electronmass\  = 9.1093818\ *10\^\(-31\);\)\)
\!\(\(permittivity\  = \ 8.854187817\ *\ 10\^\(-12\);\)\)
\!\(hbar\  = \(6.62606876*10\^\(-34\)\)\/\(2  \[Pi]\)\)
(*
\!\(1.0545715964207857`*^-34\)
*)
\!\(HartreeeV\  = \ \(electronmass\ electroncharge\^3\)\/\(16\ \(\[Pi]\^2\) permittivity\^2\ hbar\^2\)\)
(*
27.2114
*)
(*
Constants obtained from http://physics.nist.gov/cuu/Constants/ in S.I. units \
( kg m s)
*)
(*
File Loading
*)
\!\(\(\(\ \)\(Idef[0, 0, 0, 0, 0, 0, 0]\  = \ \@\(\[Pi]\/2\);\)\)\)
\!\(LoadIntegralandHamiltonianData[]\  := \ Module[{}, \[IndentingNewLine]If[Position[FileNames[\ ], "\<Ideffile\>"]\  == {}, \ Print["\<File 'Ideffile' not found in directory: \>", \ Directory[]\ , \ "\< -  File being created.\>"]; \ Idef[0, 0, 0, 0, 0, 0, 0]\  = \ \@\(\[Pi]\/2\); \ Save["\<Ideffile\>", Idef]]; \[IndentingNewLine]Print["\<The size of the Hamiltonian matrix element definition file 'Ideffile' is \>", FileByteCount["\<Ideffile\>"]\ , \ "\< Bytes.\>"]; \[IndentingNewLine]<< Ideffile;\[IndentingNewLine]\[IndentingNewLine]]\)
(*
This loads pre-calculated values for the 2 electron interaction integral \!\(\*FormBox[
  RowBox[{"\[LeftAngleBracket]", 
    RowBox[{\(\(\[Psi]\_nm\) \[Psi]\_pq\), "|", \(\[CapitalEta]\&^\_int\), "|", 
      FormBox[\(\(\[Psi]\_ab\) \[Psi]\_cd\),
        "TraditionalForm"]}], "\[RightAngleBracket]"}], TraditionalForm]\) from the file Ideffile. These values are stored as definition .e.g. \
Idef[0,0,0,0,0,0,0] = \!\(\@\(\[Pi]\/2\)\). There are only seven rather than eight quantum number inputs, because if \
the sum of the orbital angular momentum quantum numbers are not equal, i.e. m \
+q != b +d, then the interaction integral is zero. It is also a good idea to \
back-up Ideffile occasionally so lengthy calculations don't have to be \
repeated if the file is lost or corrupted. 
*)
If[Position[FileNames[ ], "alfadataexact"] == {}, 
    Print["File 'alfadataexact' not found in directory: ", Directory[], 
      " -  File being created."]; 
    Table[0, {n, 0, 2} , {m, 0, 2}, {p, 0, 2}] >> alfadataexact];

alfexact = << alfadataexact;
(*
The components of the 2-electron interaction integral (which are pieced \
together to form the values stored in Ideffile) are stored in alfadataexact \
which are then loaded into the variable "alfexact". If this file doesn't \
exist, it is created. 
*)
(*
3.3 Hamiltonian Matrix Elements
*)
(*
The functions in this section calculate the interaction Hamiltonian matrix \
using the calculations and definitions from the previous sections. 
*)
\!\(\(\(Hamileleminteract[repa_AngleBracket, \ repb_AngleBracket]\  := \ Module[{ii, j}, Which[OrbitalDifference[repa[\([1]\)], \ repb[\([1]\)]]\  >= \ 3, \ 0, \ \[IndentingNewLine]\[IndentingNewLine]\[IndentingNewLine]repa[\([1]\)]\  == \ repb[\([1]\)], \[IndentingNewLine]\ numlist\  = \ Table[ico, \ {ico, 1, NN}]; \ \[IndentingNewLine]KronOrbIdent[ii_, j_]\  := If[repb[\([1, ii]\)]\  == \ repb[\([1, j]\)], \ 1, 0]; \[IndentingNewLine]\[Sum]\+\(ii = 1\)\%\(NN - 1\)\((\[Sum]\+\(j = ii + 1\)\%NN\ \ 2\^\(-KronOrbIdent[ii, j]\)\ \((\ KroneckerDelta[repa[\([2]\)], \ repb[\([2]\)]]*\ IIIIdymexact[repa[\([1, ii, 1]\)], \ repa[\([1, ii, 2]\)], \ repa[\([1, j, 1]\)], \ repa[\([1, j, 2]\)], repb[\([1, ii, 1]\)], \ repb[\([1, ii, 2]\)], \ repb[\([1, j, 1]\)], \ repb[\([1, j, 2]\)]]\  - \ \[IndentingNewLine]\[IndentingNewLine]UUUcomp[\ P\_\(ii, j\)[numlist], \ VVV\_\(Floor[NN\/2]\), repa[\([2]\)], \ repb[\([2]\)]]\ *IIIIdymexact[repa[\([1, ii, 1]\)], \ repa[\([1, ii, 2]\)], \ repa[\([1, j, 1]\)], \ repa[\([1, j, 2]\)], repb[\([1, j, 1]\)], \ repb[\([1, j, 2]\)], \ repb[\([1, ii, 1]\)], \ repb[\([1, ii, 2]\)]])\))\), \[IndentingNewLine]\[IndentingNewLine]\[IndentingNewLine]OrbitalDifference[repa[\([1]\)], \ repb[\([1]\)]]\  == \ 2, \ \[IndentingNewLine]difforb\  = \ MultiComplement[repa[\([1]\)], \ MultiIntersection[repa[\([1]\)], repb[\([1]\)]]]; \[IndentingNewLine]If[difforb[\([1]\)]\  == \ difforb[\([2]\)], \ \[IndentingNewLine]ii\  = \ \(Position[repa[\([1]\)], \ difforb[\([1]\)]]\)[\([1, 1]\)]; \ \[IndentingNewLine]j\  = \ \(Position[repa[\([1]\)], \ difforb[\([2]\)]]\)[\([2, 1]\)], \ \[IndentingNewLine]ii\  = \ \(Position[repa[\([1]\)], \ difforb[\([1]\)]]\)[\([1, 1]\)]; \[IndentingNewLine]j\  = \ \(Position[repa[\([1]\)], \ difforb[\([2]\)]]\)[\([1, 1]\)]]; \[IndentingNewLine]LLij\  = \ LineUpPair[ii, j, repa[\([1]\)], \ repb[\([1]\)]]; \[IndentingNewLine]repc\  = \ \(repb[\([1]\)]\)[\([LLij]\)]; \[IndentingNewLine]KronOrbij\  = If[\(\(repb[\([1]\)]\)[\([LLij]\)]\)[\([ii]\)]\  == \ \(\(repb[\([1]\)]\)[\([LLij]\)]\)[\([j]\)], \ 1, 0]; \[IndentingNewLine]Nuv\  = \ \[Sqrt]2\^\(\(-\((Length[repa[\([1]\)]] - Length[Union[repa[\([1]\)]]])\)\)\(+\)\((Length[repb[\([1]\)]] - Length[Union[repb[\([1]\)]]])\)\(\ \)\); \[IndentingNewLine]numdubuv\  = \ \ If[difforb[\([1]\)]\  == \ difforb[\([2]\)], \ 0, \ \((Count[repa[\([1]\)], \ difforb[\([1]\)]] - 1)\) + \((Count[repa[\([1]\)], \ difforb[\([2]\)]] - 1)\)]; \[IndentingNewLine]\ Nuv\ 2\^\(-KronOrbij\)\ 2\^numdubuv\ \(\((\(-1\))\)\^Parity[LLij]\) \((UUUcomp[LLij\ , \ VVV\_\(Floor[NN\/2]\), repa[\([2]\)], \ repb[\([2]\)]]*\ IIIIdymexact[repa[\([1, ii, 1]\)], \ repa[\([1, ii, 2]\)], \ repa[\([1, j, 1]\)], \ repa[\([1, j, 2]\)], repc[\([ii, 1]\)], \ repc[\([ii, 2]\)], repc[\([j, 1]\)], \ repc[\([j, 2]\)]]\ \  - \ \[IndentingNewLine]\[IndentingNewLine]UUUcomp[P\_\(ii, j\)[LLij]\ , \ VVV\_\(Floor[NN\/2]\), repa[\([2]\)], \ repb[\([2]\)]]*IIIIdymexact[repa[\([1, ii, 1]\)], \ repa[\([1, ii, 2]\)], \ repa[\([1, j, 1]\)], \ repa[\([1, j, 2]\)], repc[\([j, 1]\)], \ repc[\([j, 2]\)], repc[\([ii, 1]\)], \ repc[\([ii, 2]\)]])\), \[IndentingNewLine]\[IndentingNewLine]\[IndentingNewLine]OrbitalDifference[repa[\([1]\)], \ repb[\([1]\)]]\  == \ 1, \[IndentingNewLine]difforb\  = \ MultiComplement[repa[\([1]\)], \ MultiIntersection[repa[\([1]\)], repb[\([1]\)]]]; \[IndentingNewLine]mm\  = \ \(Position[repa[\([1]\)], \ difforb[\([1]\)]]\)[\([1, 1]\)]; \[IndentingNewLine]ndu\  = \ Length[Position[repa[\([1]\)], \ difforb[\([1]\)]]]; \[IndentingNewLine]Nuv\  = \ \[Sqrt]2\^\(\(-\((Length[repa[\([1]\)]] - Length[Union[repa[\([1]\)]]])\)\)\(+\)\((Length[repb[\([1]\)]] - Length[Union[repb[\([1]\)]]])\)\(\ \)\); \[IndentingNewLine]Lmm\  = \ LineUpOne[mm, \ repa[\([1]\)], \ repb[\([1]\)]]; \[IndentingNewLine]repc\  = \ \(repb[\([1]\)]\)[\([Lmm]\)]; \[IndentingNewLine]KronOrbOne[icon2_, jcon2_]\  := If[\(\(repb[\([1]\)]\)[\([Lmm]\)]\)[\([icon2]\)]\  == \ \(\(repb[\([1]\)]\)[\([Lmm]\)]\)[\([jcon2]\)], \ 1, 0]; \[IndentingNewLine]If[repa[\([1, mm, 2]\)]\  == \ \(\(repb[\([1]\)]\)[\([Lmm]\)]\)[\([mm, 2]\)], \[IndentingNewLine]Nuv\ ndu \(\[Sum]\+\(ii = 1\)\%NN\((If[ii\  == \ mm, \ 0, \ \ \(2\^\(-KronOrbOne[ii, mm]\)\) \(\((\(-1\))\)\^Parity[Lmm]\) \((UUUcomp[Lmm, Vtrans, repa[\([2]\)], \ repb[\([2]\)]]\ *IIIIdymexact[repa[\([1, ii, 1]\)], \ repa[\([1, ii, 2]\)], \ repa[\([1, mm, 1]\)], \ repa[\([1, mm, 2]\)], repc[\([ii, 1]\)], \ repc[\([ii, 2]\)], repc[\([mm, 1]\)], \ repc[\([mm, 2]\)]]\  - \ UUUcomp[\ P\_\(ii, mm\)[Lmm], Vtrans, repa[\([2]\)], \ repb[\([2]\)]]*IIIIdymexact[repa[\([1, ii, 1]\)], \ repa[\([1, ii, 2]\)], \ repa[\([1, mm, 1]\)], \ repa[\([1, mm, 2]\)], repc[\([mm, 1]\)], \ repc[\([mm, 2]\)], repc[\([ii, 1]\)], \ repc[\([ii, 2]\)]])\)])\)\) - If[ndu\  == \ 2, \ \(\((\(-1\))\)\^Parity[Lmm]\) \((UUUcomp[Lmm, Vtrans, repa[\([2]\)], \ repb[\([2]\)]]\ *IIIIdymexact[repa[\([1, mm, 1]\)], \ repa[\([1, mm, 2]\)], \ repa[\([1, mm + 1, 1]\)], \ repa[\([1, mm + 1, 2]\)], repc[\([mm, 1]\)], \ repc[\([mm, 2]\)], repc[\([mm + 1, 1]\)], \ repc[\([mm + 1, 2]\)]]\  - \ UUUcomp[\ P\_\(mm, mm + 1\)[Lmm], Vtrans, repa[\([2]\)], \ repb[\([2]\)]]*IIIIdymexact[repa[\([1, mm, 1]\)], \ repa[\([1, mm, 2]\)], \ repa[\([1, mm + 1, 1]\)], \ repa[\([1, mm + 1, 2]\)], \(\(repb[\([1]\)]\)[\([\ Lmm]\)]\)[\([mm + 1, 1]\)], \ \(\(repb[\([1]\)]\)[\([\ Lmm]\)]\)[\([mm + 1, 2]\)], \(\(repb[\([1]\)]\)[\([\ Lmm]\)]\)[\([mm, 1]\)], \ \(\(repb[\([1]\)]\)[\([\ Lmm]\)]\)[\([mm, 2]\)]])\)\ \ Nuv, 0], 0]\[IndentingNewLine]]]\)\(\[IndentingNewLine]\)
  \)\)
(*
\!\(\*
  RowBox[{\(General::"spell"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(repc\\)\\\" is similar to existing symbols \\!\\({repa, repb}\\). \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell\\\"]\\)\"\>"}]\)
*)
(*
Hamileleminteract is a implementation of the rules for the calculation of the \
interaction matrix element between two basis functions.  Hamileleminteract \
breaks down cases on the basis of the "orbital difference" of the two basis \
elements and then applies the formula appropriate to the case. 
*)
\!\(HamiltonianSparseLsortInt[basisgen_List]\  := \ Module[{tempelem}, \ \[IndentingNewLine]basissort\  = \ Split[Sort[basisgen, \((\ \[Sum]\+\(i = 1\)\%NN\((#1[\([1, i, 2]\)])\)\  < \ \ \[Sum]\+\(i = 1\)\%NN\((#2[\([1, i, 2]\)])\))\) &], \ \[Sum]\+\(i = 1\)\%NN\((#1[\([1, i, 2]\)])\)\  == \ \ \[Sum]\+\(i = 1\)\%NN\((#2[\([1, i, 2]\)])\) &]; \[IndentingNewLine]numL\  = \ Length[basissort]; \[IndentingNewLine]lenhist\  = \ Length\  /@ \ basissort; \[IndentingNewLine]Hlist\  = {}; \[IndentingNewLine]For[kco\  = 1, \ kco <= \ numL, \[IndentingNewLine]hist\  = \ Plus\  @@ \ lenhist[\([Range[1, kco - 1]]\)]; \[IndentingNewLine]For[ico = 1, \ ico <= lenhist[\([kco]\)], \ \[IndentingNewLine]For[jco = ico, jco <= lenhist[\([kco]\)], \ \ \[IndentingNewLine]tempelem\  = \ Hamileleminteract[basissort[\([kco, ico]\)], \ basissort[\([kco, jco]\)]]; \[IndentingNewLine]If[tempelem\  == \ 0, , \ Hlist\  = \ Join[Hlist, \ {{hist\  + \ ico, hist\  + \ jco} -> \ tempelem, \ {hist + jco, hist\  + \ ico} -> \ tempelem}]]; \(jco++\)]; \[IndentingNewLine]\(ico++\)]; \ \[IndentingNewLine]\(kco++\)]\ ; SparseArray[Hlist]]\)
\!\(Triag[x_]\  := \ \(1\/2\) x \((x + 1)\)\)
(*
\!\(\*
  RowBox[{\(General::"spell"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(Triag\\)\\\" is similar to existing symbols \\!\\({Trial, Trig}\\). \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell\\\"]\\)\"\>"}]\)
*)
\!\(Progress[ico2_, \ len2_]\  := \ ico2\/2 - ico2\^2\/2 + ico2\ len2\)
HamiltonianInt[basisgen_List] := Module[{tempelem, H, tot}, 
    len = Length[basisgen];
    tot = 0;
    H = SparseArray[{}, {len, len}];
    For[ico = 1, ico <= len, 
      If[Floor[Progress[ico, len]/Triag[len] *20] > tot, 
        tot = Floor[Progress[ico, len]/Triag[len] *20];  Print[tot*5, "%"]]; 
      For[jco = ico, jco <= len, 
        H [[ico, jco]] = 
          H[[jco, ico]] = Hamileleminteract[basisgen[[ico]], basisgen[[jco]]];
         jco++]; ico++] ; H]
(*
HamiltonianInt and HamiltonianSparseLsortInt calculate the full interaction \
Hamiltonian matrix for a given basis. HamiltonianInt is the most commonly \
used but HamiltonianSparseLsortInt is slightly faster when the basis used \
contains elements of different values of the total orbital angular momentum. \
The calculation of the interaction Hamiltonian matrix is usually the most \
time intensive step.
*)
HamiltonianZeroMagnetic[basisgen_List, firengap_, magneticfield_, relmass_, 
    relperm_] := 
  SparseArray[
    Table[{i, i} -> 
        NonInteractingEnergyMagnetic[basisgen[[i, 1]], firengap, 
          magneticfield, relmass, relperm], {i, 1, 
        Length[basisgen]}], {Length[basisgen], Length[basisgen]}]
(*
HamiltonianZeroMagnetic gives the non-interacting Hamiltonian matrix for \
specific values of physical parameters.  This matrix can then be combined \
with the interaction matrix to give the final Hamiltonian matrix from which \
energy levels can be calculated.
*)
\!\(Hfinal[Hint_SparseArray, \ basisgen_List, \ firengap_, \ magfield_, \ relmass_, \ relperm_]\  := \ \@kfunction[firengap, magfield, \ relmass\ , \ relperm]\ Hint\  + \ HamiltonianZeroMagnetic[\ basisgen, \ firengap, magfield, \ relmass, \ relperm]\)
\!\(HfinalSorted[Hint_SparseArray, \ basisgen2_List, \ firengap_, \ magfield_, \ relmass_, \ relperm_]\  := \ \@kfunction[firengap, magfield, \ relmass, \ relperm]\ Hint\  + \ HamiltonianZeroMagnetic[\ Sort[basisgen2, \((\ \[Sum]\+\(i = 1\)\%NN\((#1[\([1, i, 2]\)])\)\  < \ \ \[Sum]\+\(i = 1\)\%NN\((#2[\([1, i, 2]\)])\))\) &], \ firengap, magfield, \ relmass, \ relperm]\)
(*
Hfinal and HfinalSorted give the final Hamiltonian matrix from which energy \
levels can be calculated.  The former is for use for interaction matrices \
calculated using HamiltonianInt while the latter is for use for interaction \
matrices calculated using HamiltonianSparseLsortInt.
*)
InteractionHamiltonianMatrixElement[repa_AngleBracket, repb_AngleBracket] := 
  Hamileleminteract[repa, repb]
InteractionHamiltonian[basisgen_List] :=  HamiltonianInt[basisgen]
TotalHamiltonian[Hint_SparseArray, basisgen_List, firengap_, magfield_, 
    relmass_, relperm_] := 
  Hfinal[Hint, basisgen, firengap, magfield, relmass, relperm] 
(*
3.4 Two Electron Interaction Function Definitions
*)
(*
Intermediate integral
*)
\!\(alphadiagexact[\[Alpha]_, \[Beta]_]\  := \ 2\^\(\(-3\)\ \[Alpha] - \[Beta]\)\ \[Pi]\^\(3/2\)\ Binomial[2\ \[Alpha], \[Alpha]]\ Gamma[1 + \[Alpha]]\ Gamma[1\/2\ \((1 + \[Alpha] + \[Beta])\)]\ HypergeometricPFQRegularized[{1\/4\ \((1 + 2\ \[Alpha])\), 1\/4\ \((3 + 2\ \[Alpha])\), 1\/2\ \((1 + \[Alpha] + \[Beta])\)}, {1 + \[Alpha], 1\/2\ \((2 + \[Alpha] + \[Beta])\)}, 1]\)
\!\(FF\_0[\[Beta]_]\  = \ eff[\[Beta]]\)
(*
eff[\[Beta]]
*)
\!\(FF\_1[\[Beta]_]\  = \ eff[\[Beta]]\/2\)
(*
\!\(eff[\[Beta]]\/2\)
*)
\!\(FF\_2[\[Beta]_]\  = \ eff[\[Beta]]\/2 - \ eff[\[Beta] + 2]\)
(*
\!\(eff[\[Beta]]\/2 - eff[2 + \[Beta]]\)
*)
\!\(FF\_n_[\[Beta]_]\  := \ \(FF\_n[\[Beta]]\  = \ \(\(FF\_\(n - 1\)[\[Beta]]\)\(\ \)\(-\)\(FF\_\(n - 2\)[\[Beta] + 2]\)\(\ \)\)\)\)
\!\(alfacalcexact[\[Alpha]_, \ a_, \ b_]\  := \ FF\_\(Abs[b - a]\/2\)[Min[a, b]]\  /. \ eff[x_]\  -> \ alphadiagexact[\[Alpha] - 1, x]\)
(*
The alfacalcexact functions calculate the value of the integral\!\(TraditionalForm\`\[Integral]\_0\%\(\[Pi]/2\)\(\[Integral]\_0\%\(2  \[Pi]\)\(\(Cos[\[Alpha]\[Delta]] \(Sin[\[Lambda]]\^a\) Cos[\[Lambda]]\^b\)\/\@\(1 - Sin[2  \[Lambda]] Cos[2  \[Lambda]]\)\) \[DifferentialD]\[Delta]\)\)\[DifferentialD]\[Lambda] for different values of \[Alpha], a and b. These \
arise in the calculation of the 2-electron interaction integral. 
*)
alfadymexact[\[Alpha]_, a_, b_] := Module[{dim = Dimensions[alfexact]},
    If[\[Alpha] > Dimensions[alfexact][[1]], 
      alfexact = 
        PadRight[
          alfexact, {\[Alpha], Dimensions[alfexact][[2]], 
            Dimensions[alfexact][[3]]}], Null ];
    If[a > Dimensions[alfexact][[2]], 
      alfexact = 
        PadRight[
          alfexact, {Dimensions[alfexact][[1]], a, 
            Dimensions[alfexact][[3]]}], Null ];
    If[b > Dimensions[alfexact][[3]], 
      alfexact = 
        PadRight[
          alfexact, {Dimensions[alfexact][[1]], Dimensions[alfexact][[2]], 
            b}], Null  ];
    If[alfexact[[\[Alpha], a, b]] == 0, 
      alfexact[[\[Alpha], a, b]] = alfacalcexact[\[Alpha], a, b], 
      alfexact[[\[Alpha], a, b]]]]
(*

*)
(*
The cell above calculates, stores and retrieves portions of the two electron \
interaction which are later pieced together to give the whole integral.
*)
(*
Interaction Integral
*)
\!\(LLL\_\(n_, m_\)[r_]\  := \ Binomial[n + Abs[m], n]\ Hypergeometric1F1[\(-n\), Abs[m] + 1, \ r]\)
\!\(R\_\(n_, m_\)[r_, \ k_] := \ \((2\ \(k\^\(Abs[m] + 1\)\) \(n!\)\/Gamma[Abs[m] + n + 1])\)\^\(1\/2\)\ \(r\^Abs[m]\) \(\[ExponentialE]\^\(\(\(-k\)\/2\) r\^2\)\) LLL\_\(n, m\)[k\ r\^2]\)
\!\(RRR\_\(n_, m_\)[r_] := \ \((2\ \(n!\)\/Gamma[Abs[m] + n + 1])\)\^\(1\/2\)\ \(r\^Abs[m]\) LLL\_\(n, Abs[m]\)[r\^2]\)
(*
LLL gives the generalized Laguerre polynomials, R gives the Radial \
wavefunction, and RRR is the radial wave function with the exponential term \
removed. 
*)
\!\(GGGG[r_, s_, n_, m_, \ p_, q_, \ a_, \ b_, \ c_, \ d_]\  := \ If[Head[Expand[r\ s \( 1\/\(4  \[Pi]\^2\)\) RRR\_\(n, m\)[r]\ RRR\_\(p, q\)[s]\ RRR\_\(a, b\)[r]\ RRR\_\(c, d\)[s]]]\  === \ Plus, \[IndentingNewLine]ReplacePart[Expand[r\ s \( 1\/\(4  \[Pi]\^2\)\) RRR\_\(n, m\)[r]\ RRR\_\(p, q\)[s]\ RRR\_\(a, b\)[r]\ RRR\_\(c, d\)[s]], \ List, 0], {r\ s \( 1\/\(4  \[Pi]\^2\)\) RRR\_\(n, m\)[r]\ RRR\_\(p, q\)[s]\ RRR\_\(a, b\)[r]\ RRR\_\(c, d\)[s]}]\)
\!\(JJJJ[n_, m_, \ p_, q_, \ a_, \ b_, \ c_, \ d_]\  := \((GGGG[s, r, n, m, p, q, a, b, c, d])\) /. \ {x_\ \(r\^y_\) s\^z_\  -> \ If[y <= \ z, \ {x, y, z}, {x, z, y}], \ x_\ r\ s\^z_\  -> \ {x, 1, z}, \ x_\ r\^y_\ s\  -> \ {x, 1, y}, \ x_\ r\ s\  -> \ {x, 1, 1}}\)
\!\(IIIIcalcexact[n_, m_, \ p_, q_, \ a_, \ b_, \ c_, \ d_]\  := \ Module[{JJ = \ JJJJ[n, m, p, q, a, b, c, d]}, \[Sum]\+\(i = 1\)\%\(Length[JJ]\)JJ[\([i, 1]\)]\ alfadymexact[Abs[b - m] + 1, JJ[\([i, 2]\)], JJ[\([i, 3]\)]]\ \[Pi]\ Gamma[\(1\/2\) \((1 + JJ[\([i, 2]\)] + JJ[\([i, 3]\)])\)]]\)
IIIIdymexact[n_Integer, m_Integer, p_Integer, q_Integer, a_Integer, b_Integer,
     c_Integer, d_Integer] := 
  If[b + d == m + q, 
    If[ValueQ[Idef[n, m, p, q, a, b, c]], Idef[n, m, p, q, a, b, c], 
      Idef[n, m, p, q, a, b, c] = 
        Idef[a, b, c, d, n, m, p] = 
          Idef[p, q, n, m, c, d, a] = 
            Idef[c, d, a, b, p, q, n] = 
              IIIIcalcexact[n, m, p, q, a, b, c, d] ], 0]
(*
IIII functions calculate, store and retrieve the integral \!\(\*FormBox[
  RowBox[{"\[LeftAngleBracket]", 
    RowBox[{\(\(\[Psi]\_nm\) \[Psi]\_pq\), "|", \(\[CapitalEta]\&^\_int\), "|", 
      FormBox[\(\(\[Psi]\_ab\) \[Psi]\_cd\),
        "TraditionalForm"]}], "\[RightAngleBracket]"}], TraditionalForm]\). 
*)
InteractionIntegral[n_Integer, m_Integer, p_Integer, q_Integer, a_Integer, 
    b_Integer, c_Integer, d_Integer] := IIIIdymexact[n, m, p, q, a, b, c, d] 
(*
3.5 Additional Functions
*)
(*
Conversion Factors
*)
\!\(energyconversionfactor[relmass_, \ relperm_]\  := \ \(relmass\/relperm\^2\) HartreeeV\)
\!\(lengthscale[relmass_, \ relperm_]\  := \ \(4  \[Pi]\ relperm\ permittivity\ hbar\^2\)\/\(relmass\ electronmass\ electroncharge\^2\)\ *10\^9\)
(*
The numerical computations are done in "material units", which are atomic \
units adjusted to the permittivity and effective electron mass of the \
material we are working with. To convert to ordinary units, energy values are \
multiplied by energyconversionfactor and length values are multiplied by \
lengthscale to convert to eV and nm respectively. Note that lengthscale is \
the effective Bohr radius.
*)
\!\(t\_a[relmass_, \ relperm_]\  := \ hbar/\((relmass*\ electronmass*\((lengthscale[relmass, \ relperm]*10\^\(-9\))\)\^2\ )\)\)
\!\(Bconversion[relmass_, \ relperm_] := electroncharge\/\(relmass*\ electronmass*\ t\_a[relmass, \ relperm]\)\)
(*
Bconversion is the conversion factor between the magnetic field measured in \
Tesla to material units. i.e. 1 T = Bconversion in material units. Similarly, \!\(\*FormBox[
  StyleBox[\(t\_\(\(a\)\(\ \)\)\),
    FontWeight->"Bold"], TraditionalForm]\) is the number of material time units in one second.
*)
\!\(\(\[Omega]0[firengap_, \ relmass_, \ relperm_]\  := \ firengap\/\(\(1000\)\(\ \)\(HartreeeV\)\(\ \)\)\ *\ relperm\^2\/\(\(\ \)\(relmass\)\);\)\)
(*
\[Omega]0 gives the value of the harmonic well constant that gives the \
specified energy level spacing 
*)
\[Omega]c[magfield_, relmass_, relperm_] := 
  magfield*Bconversion[relmass, relperm]
\!\(lnought[magfield_, \ relmass_, \ relperm_] := \@\(\(\ \)\(1\/\(magfield*Bconversion[relmass, \ relperm]\)\)\)\)
\!\(bfunction[firengap_, \ magfield_, \ relmass_, \ relperm_]\  := \ \@\(1\  + 4 \((\ \[Omega]0[firengap, \ relmass, \ relperm]\/\[Omega]c[magfield, \ relmass, \ relperm])\)\^2\)\)
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(bfunction\\)\\\" is similar to existing symbol \\\"\\!\\(kfunction\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
\!\(kfunction[firengap_, \ magfield_, \ relmass_, \ relperm_] := \ If[magfield\  == \ 0, \[Omega]0[firengap, \ relmass, \ relperm], \ bfunction[firengap, \ magfield, \ relmass, \ relperm]\/\(2\ lnought[magfield, \ relmass, \ relperm]\^2\)]\)
(*
The above cells are used to handle the conversion between atomic/material \
units in which the calculations are performed and units convenient for input \
and output.
*)
(*
List Manipulation
*)
MultiIntersection[l1_List, l2_List] := 
  Module[{nl, f}, f[x_] := {First[#], Length[#]} & /@ Split[Sort[x]]; 
    nl = Sort[Join[Flatten[Map[f, {l1, l2}], 1]]]; 
    nl = Split[nl, #[[1]] === #2[[1]] &]; 
    Flatten[Cases[nl, {{x_, m_}, {x_, n_}} :> Table[x, {m}]], 1]]
(*
MultiIntersection gives the intersection of two lists with any element \
appearing more than once in both lists also appearing more than once in the \
intersection (corresponding to the minimum number of times it appears in \
either list).
*)
MultiComplement[l1_List, l2_List] := Module[{nl, f},
    f[x_] := {First[#], Length[#]} & /@ Split[Sort[x]];
    nl = Sort[Join[Flatten[Map[f, {l1, l2}], 1]]];
    nl = Split[nl, #[[1]] === #2[[1]] &];
    Flatten[Cases[nl, {{x_, m_}, {x_, n_}} :> Table[x, {Abs[n - m]}]], 1]~
      Join~Flatten[Cases[nl, {{x_, m_}} :> Table[x, {Abs[m]}]], 1]]
(*
MultiComplement gives the elements of the first list that don't appear in the \
second a corresponding number of times.
*)
OrbitalDifference[u_List, v_List] := 
  Length[u] - Length[MultiIntersection[u, v]]
(*
OrbitalDifference gives the number of orbitals by which two spatial \
wavefunctions differ. If this is greater than or equal to three then the \
interaction between the two spatial wavefunctions is zero.
*)
LineUpPair[ilineup_, jlineup_, u_List, v_List] := Module[{postemp},
    dubpos = Table[0, {ico, 1, Length[u]}];
    If[ OrbitalDifference[Delete[u, {{ilineup}, {jlineup}}], v] == 0,
      
      tempaline = Table[
          
          Which[
            k == ilineup, ihold, 
            k == jlineup, jjhold,
             1 == 1, postemp = Flatten[Position[v, u[[k]]]]; 
            
             If[Length[postemp] == 1, 
              postemp[[1]], 
              
               dubpos[[ Position[u, u[[k]]] [[1, 1]] ]]++; 
              postemp[[  dubpos[[Position[u, u[[k]]][[1, 1]]]]] ]]
            ], {k, 1, Length[u]}];  
      ijreplace = 
        Complement[Table[ico, {ico, 1, Length[u]}], 
          Delete[tempaline, {{ilineup}, {jlineup}}]]; 
      tempaline /. {ihold -> ijreplace[[1]], jjhold -> ijreplace[[2]]},
      
      0] ]
(*
LineUpPair gives the permutation that line up all the orbitals in u and v \
except the orbitals in positions i and j.  If it can't find one then it \
returns 0.
*)
LineUpOne[jlineup_, u_List, v_List] := Module[{postemp},
    dubpos = Table[0, {ico, 1, Length[u]}];
    If[ OrbitalDifference[Delete[u, {{jlineup}}], v] == 0,
      tempaline = Table[
          If[ k == jlineup, jjhold, 
            postemp = Flatten[Position[v, u[[k]]]]; 
            If[Length[postemp] == 1, 
              postemp[[1]], 
               dubpos[[ Position[u, u[[k]]] [[1, 1]] ]]++; 
              postemp[[  dubpos[[Position[u, u[[k]]][[1, 1]]]]] ]]],
          {k, 1, Length[u]}]; 
       ijreplace = 
        Complement[Table[ico, {ico, 1, Length[u]}], 
          Delete[tempaline, jlineup]]; tempaline /. jjhold -> ijreplace[[1]], 
      0] ]
(*
LineUpOne gives the permutation that line up all the orbitals in u and v \
except the orbital in position j. If it can't find one then it returns 0.
*)
(*
Energy and Parity 
*)
\!\(SingleOrbitalEnergyMagnetic[u_List, \ firengap_, \ magfield_, \ relmass_, \ relperm_]\  := \ \((2\ u[\([1]\)]\  + \ Abs[u[\([2]\)]]\  + 1)\)*\ \[Sqrt]\((\((\[Omega]0[firengap, \ relmass, \ relperm])\)\^2 + \ \(1\/4\) \((\[Omega]c[magfield, \ relmass, \ relperm])\)\^2)\)\ \  - \ \(u[\([2]\)]\/2\) \[Omega]c[magfield, \ relmass, \ relperm]\)
(*
SingleOrbitalEnergyMagnetic gives the energy eigenvalue of a single electron \
in a given harmonic well under the influence of a given magnetic field, where \
the spatial wavefunction is specified by the quantum numbers in u.
*)
NonInteractingEnergyMagnetic[uli_List, firengap_, magfield_, relmass_ , 
    relperm_] := 
  Plus @@ (SingleOrbitalEnergyMagnetic[#1, firengap, magfield, relmass, 
            relperm] & /@ uli)
(*
NonInteractingEnergyMagnetic gives the energy of multiple non interaction \
electrons in a given harmonic well under the influence of a given magnetic \
field with the spatial wavefunctions  specified by the quantum numbers in \
uli.
*)
SingleOrbitalEnergy[v_List] := 2v[[1]] + Abs[v[[2]]] + 1
(*
SingleOrbitalEnergy gives the energy eigenvalue of a single electron in a \
given harmonic well at zero magnetic field where the spatial wavefunction is \
specified by the quantum numbers in v.
*)
NonInteractingEnergy[uli_List] := Plus @@ SingleOrbitalEnergy /@ ( uli)
(*
NonInteractingEnergy gives the energy of multiple non interaction electrons \
in a given harmonic well at zero magnetic field with the spatial \
wavefunctions again specified by the quantum numbers in uli.
*)
Parity[u_List] := 
  Mod[PermutationType[u] . Table[Mod[i + 1, 2], {i, 1, Length[u]}], 2]
(*
Parity gives the Parity of a permutation, 1 for odd and 0 for even.
*)
(*
3.6 Wavefunctions and Visualization
*)
\!\(\(SingleElectronProbablityDensity[n_, m_, \ firstenergygap_, \ magneticfield_, \ relativeeffectivemass_, \ relativepermittivity_, \ rmax_, opts___]\  := \ ParametricPlot3D[{r\ lengthscale[relativeeffectivemass, \ relativepermittivity] Cos[\[Alpha]], \ r\ lengthscale[relativeeffectivemass, \ relativepermittivity]\ Sin[\[Alpha]], \ \((R\_\(n, m\)[r, \ kfunction[firstenergygap, \ magneticfield, \ relativeeffectivemass, \ relativepermittivity]])\)\^2}, \ {r, 0.0001, rmax}, \ {\[Alpha], 0, 2  \[Pi]}, \ opts];\)\)
\!\(SingleExtractor[repa_AngleBracket, \ repb_AngleBracket, var_, angle_, \ k_]\  := \ Module[{}, \ \[IndentingNewLine]Which[OrbitalDifference[repa[\([1]\)], repb[\([1]\)]] >= \ 2, 0, \[IndentingNewLine]\[IndentingNewLine]OrbitalDifference[repa[\([1]\)], repb[\([1]\)]] == 1, \[IndentingNewLine]\[IndentingNewLine]difforb\  = \ MultiComplement[repa[\([1]\)], \ MultiIntersection[repa[\([1]\)], repb[\([1]\)]]]; \[IndentingNewLine]mm\  = \ \(Position[repa[\([1]\)], \ difforb[\([1]\)]]\)[\([1, 1]\)]; \[IndentingNewLine]difforbtwo\  = \ MultiComplement[repb[\([1]\)], \ MultiIntersection[repa[\([1]\)], repb[\([1]\)]]]; \[IndentingNewLine]pp\  = \ \(Position[repb[\([1]\)], \ difforbtwo[\([1]\)]]\)[\([1, 1]\)]; \[IndentingNewLine]Nuv\  = \ \[Sqrt]2\^\(\(-\((Length[repa[\([1]\)]] - Length[Union[repa[\([1]\)]]])\)\)\(+\)\((Length[repb[\([1]\)]] - Length[Union[repb[\([1]\)]]])\)\(\ \)\); \[IndentingNewLine]Lmm\  = \ LineUpOne[mm, \ repa[\([1]\)], \ repb[\([1]\)]]; \[IndentingNewLine]\((\(Nuv\/NN\) \(\((\(-1\))\)\^Parity[Lmm]\) UUUcomp[Lmm\ , \ VVV\_\(Floor[NN\/2]\), repa[\([2]\)], \ repb[\([2]\)]]*\(\(1\/\(\(2\)  \(\[Pi]\)\(\ \)\)\) \((R\_\(difforb[\([1, 1]\)], difforb[\([1, 2]\)]\)[var, \ k]\ R\_\(difforbtwo[\([1, 1]\)], difforbtwo[\([1, 2]\)]\)[var, \ k])\)\) \[ExponentialE]\^\(\[ImaginaryI] \((difforbtwo[\([1, 2]\)] - difforb[\([1, 2]\)])\)\ angle\))\), \[IndentingNewLine]OrbitalDifference[repa[\([1]\)], repb[\([1]\)]] == 0, \[IndentingNewLine]\[IndentingNewLine]\[Sum]\+\(orbcount = 1\)\%NN\((R\_\(repa[\([1, orbcount, 1]\)], repb[\([1, orbcount, 2]\)]\)[var, \ k]\ )\)^2\/\(\(2\)  \(\[Pi]\)\(\ \)\)]]\)
\!\(ElectronDensity[vec_List, \ basis_List, var_, angle_, \ k_]\  := \ \[Sum]\+\(icoun = 1\)\%\(Length[vec]\)\(\[Sum]\+\(jcoun = 1\)\%\(Length[vec]\)If[vec[\([icoun]\)] == 0\  || \ vec[\([jcoun]\)] == 0, \ 0, vec[\([icoun]\)] vec[\([jcoun]\)]\ SingleExtractor[basis[\([icoun]\)], \ basis[\([jcoun]\)], var, angle, \ k]]\)\)
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(icoun\\)\\\" is similar to existing symbol \\\"\\!\\(icount\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(jcoun\\)\\\" is similar to existing symbol \\\"\\!\\(icoun\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
(*
ElectronDensity allows the calculation of the electron density function from \
eigenvectors calculated from the total Hamiltonian matrix.
*)
EnergyLevelPlot[Hfin_SparseArray, numberoflevels_, relmass_, relperm_,  
    opts___] := 
  ListPlot[Reverse[
        Eigenvalues[
            Hfin, -numberoflevels][[Range[1, 
              numberoflevels]]]] energyconversionfactor[relmass, relperm], 
    opts]
(*
EnergyLevelPlot plots the energy levels for the "n" lowest levels where "n" \
is defined by the input  numberoflevels.
*)
Convergence[Hint_, basis_, firengap_, magfield_, relmass_, relperm_, Lvalue_, 
    sizelist_List, levlist_] := Module[{H, enlev = {}, basissizelist = {}},
    For[aco = 1, aco <= Length[sizelist],
      H = 
        Hfinal[Hint[[Range[1, sizelist[[aco]]], Range[1, sizelist[[aco]]]]], 
          basis[[Range[1, sizelist[[aco]]]]], firengap, magfield, relmass, 
          relperm];
      enlev = 
        Append[enlev, 
          energyconversionfactor[relmass, 
              relperm] Reverse[Eigenvalues[H, -Max[levlist]]][[levlist]]]; 
      aco++];
    enlev]

(*
Convergence produces a table of the energy levels defined by the input \
levlist versus various basis sizes determined by sizelist. This can later be \
plotted.
*)
LevelVsField[Hint_, basis_, firengap_, relmass_, relperm_, Lvalue_ , magmin_, 
    magmax_, step_ , levlist_] := Module[{H, enlev = {}},
    For[bco = magmin, bco <= magmax,
      H = Hfinal[Hint, basis, firengap, bco, relmass, relperm];
      enlev = 
        Append[enlev, {bco}~
            Join~ (energyconversionfactor[relmass, 
                  relperm]Reverse[Eigenvalues[H, -Max[levlist]]][[levlist]])];
       bco = bco + step];
    enlev]
(*
LevelVsField produces table of the variation of the energy levels defined by \
levlist versus changes in the magnetic field.  
*)
ElectronDensityFunction[vec_, basisgen_, r_, \[Alpha]_, firstenergygap_, 
    magneticfield_, relativeeffectivemass_, relativepermittivity_] := 
  ElectronDensity[vec, basisgen, r, \[Alpha], 
    kfunction[firstenergygap, magneticfield, relativeeffectivemass, 
      relativepermittivity]]
(*
3.7 Spin Matrix Calculation
*)
(*
This section calculates spin eigenfunctions and the corresponding \
representation matrices that are used in subsequent calculations.
*)
\!\(P\_\(i_, j_\)[u_List]\  := \ Delete[Insert[Delete[Insert[u, \ u[\([j]\)], i], i + 1], \ u[\([i]\)], \ j], \ j + 1]\)
(*
\!\(\*FormBox[
  RowBox[{
    StyleBox[\(P\_\(i, j\)\),
      FontWeight->"Bold"], " "}], TraditionalForm]\) swaps the \!\(TraditionalForm\`i\^th\)and \!\(TraditionalForm\`j\^\(\(th\)\(\ \)\)\)entries in a list. Each basic eigenfunction of \!\(TraditionalForm\`S\_z\) can be represented uniquely by converting it to a decimal number by \
considering the 1's and 0's to form a binary number.  Ord is a function that \
takes a given decimal number calculated from an eigenfunction of \!\(TraditionalForm\`\(\(S\_z\)\(\ \)\)\)and gives its position in the list \[Theta].  This is a convenient way of \
retrieving the eigenfunctions when required.
*)
PPP[n_] := Module[{}, UU = ZeroMatrix[NN]; 
    For[i = 1, i <= NN, UU[[i, PP[[n, i]]]] = 1; i++]; UU]
 UUU[m_, VV_List] := Module[{f, a, Z}, 
    f = Length[VV]; 
    a = {};
    For[n = 1, n <= PermutationNumber, 
      a = Append[a, Ord[FromDigits[PPP[m].\[Theta][[n]], 2]]]; n++]; 
    UU = ZeroMatrix[PermutationNumber];
    For[i = 1, i <= PermutationNumber, UU[[i, a[[i]]]] = 1; i++];
    Z = ZeroMatrix[f];
    For[k = 1, k <= f,
      For[j = 1, j <= f, Z[[k, j]] = VV[[j]].(UU. VV[[k]]); j++]; k++]; 
    Simplify[Z]]
UUUcomp[perm_List, VV_List, irow_, jcol_] := 
  Module[{f = Length[VV], a = {}}, 
    f = Length[VV]; 
    a = {};
    For[n = 1, n <= PermutationNumber, 
      a = Append[a, Ord[FromDigits[\[Theta][[n]][[perm]], 2]]]; n++]; 
    VV[[jcol]].(VV[[irow]][[a]])]
(*
PPP[n] gives the permutation matrix corresponding to the \!\(TraditionalForm\`n\^\(\(th\)\(\ \)\)\)permutation in the list PP.
UUU gives the representation matrix of  a particular permutation.
*)
\!\(\(\(Trans[n_]\)\(\ \)\(:=\)\(\ \)\(\(Position[PP, \ P\_\(2  n - 1, 2  n\)[Table[x, {x, 1, NN}]]]\)[\([1, 1]\)]\)\(\ \)\)\)
\!\(CalculateSpinBasis[]\  := \ Module[{}, \[IndentingNewLine]NN = \ Length[ElectronList]; \[IndentingNewLine]\[Theta] = Permutations[ElectronList]; \[IndentingNewLine]PermutationNumber\  = \ Length[\[Theta]]; \[IndentingNewLine]aint\  = \ {}; \[IndentingNewLine]For[i = 1, i <= \ PermutationNumber, aint = \ Append[aint, {FromDigits[\(Permutations[ElectronList]\)[\([i]\)], \ 2], i}]; \ \(i++\)]; \ \[IndentingNewLine]Ord\  = \ Interpolation[aint, \ InterpolationOrder -> \ Min[3, NN]]; \[IndentingNewLine]MM\  = \ \(NN\/4\) \((4 - NN)\)\ IdentityMatrix[PermutationNumber]; \[IndentingNewLine]For[k = 1, \ k <= \ PermutationNumber, \ For[j = 2, j <= NN, For[i = 1, i < j, \ MM[\([Ord[FromDigits[P\_\(i, j\)[\[Theta][\([k]\)]], 2]], k]\)]\  = MM[\([Ord[FromDigits[P\_\(i, j\)[\[Theta][\([k]\)]], 2]], \ k]\)]\  + \ 1; \(i++\)]\ ; \(j++\)]; \ \(k++\)]; \[IndentingNewLine]Vectors1\  = \ NullSpace[MM - S \((S + 1)\) IdentityMatrix[PermutationNumber]]; \[IndentingNewLine]MMSparse = \ SparseArray[{i_, i_} -> \ \(NN\/4\) \((4 - NN)\) - S \((S + 1)\), {PermutationNumber, PermutationNumber}]; \[IndentingNewLine]For[k = 1, \ k <= \ PermutationNumber, \ For[j = 2, j <= NN, For[i = 1, i < j, \ MMSparse[\([Ord[FromDigits[P\_\(i, j\)[\[Theta][\([k]\)]], 2]], k]\)]\  = MMSparse[\([Ord[FromDigits[P\_\(i, j\)[\[Theta][\([k]\)]], 2]], \ k]\)]\  + \ 1; \(i++\)]\ ; \(j++\)]; \ \(k++\)]; \[IndentingNewLine]Vectors1\  = \ NullSpace[MMSparse]; \[IndentingNewLine]VV\  = \ GramSchmidt[Vectors1]; \[IndentingNewLine]PP\  = \ MinimumChangePermutations[Table[x, {x, 1, NN}]]; \[IndentingNewLine]\[IndentingNewLine]\[IndentingNewLine]Householder[{}]\  = {}; \[IndentingNewLine]numpairs\  = Floor[NN\/2]; \[IndentingNewLine]VVV\_n_\  =. ; \[IndentingNewLine]For[i = 0, \ i <= \ Floor[NN\/2], VVV\_i =. \ ; \(i++\)]\ ; \[IndentingNewLine]VVV\_0 = \ VV; \[IndentingNewLine]VVV\_\(\(n_\)\(\ \)\) := \ \(VVV\_\(\(n\)\(\ \)\) = \ Simplify[\((Householder[NullSpace[UUU[Trans[n], VVV\_\(n - 1\)] + IdentityMatrix[Length[VVV\_\(n - 1\)]]]]~Join~Householder[\ NullSpace[UUU[Trans[n], VVV\_\(n - 1\)] - IdentityMatrix[Length[VVV\_\(n - 1\)]]]])\) . VVV\_\(n - 1\)]\); \[IndentingNewLine]For[i\  = 1, \ i <= \ Floor[NN\/2], \ VVV\_i; \(i++\)]; \[IndentingNewLine]Vtrans\ \  = \ VVV\_\(Floor[NN\/2]\)\ ; \[IndentingNewLine]Norsym\  = \ Module[{imat, jmat, \ nortable, Len\  = \ Length[Vtrans]}, \[IndentingNewLine]nortable\  = \ Table[0, {numpairs}, {Len}]; \[IndentingNewLine]For[jmat = 1, \ jmat <= \ Len, \[IndentingNewLine]nortable[\([1, jmat]\)]\  = \ 1\/2 - \(1\/2\) UUUcomp[P\_\(1, 2\)[Table[x, {x, 1, NN}]], \ Vtrans, jmat, jmat\ ]; \[IndentingNewLine]\(jmat++\)]; \[IndentingNewLine]For[imat = 2, \ imat\  <= \ numpairs, \[IndentingNewLine]For[jmat = 1, \ jmat <= \ Len, \[IndentingNewLine]nortable[\([imat, jmat]\)]\  = \ nortable[\([imat - 1, jmat]\)] \((1\/2 - \(1\/2\) UUUcomp[P\_\(2  imat - 1, 2  imat\)[Table[x, {x, 1, NN}]], \ Vtrans, jmat, jmat\ ])\); \[IndentingNewLine]\(jmat++\)]; \[IndentingNewLine]\(imat++\)]; nortable];\[IndentingNewLine]]\)
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(jmat\\)\\\" is similar to existing symbol \\\"\\!\\(imat\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
(*
\!\(\*
  RowBox[{\(General::"spell1"\), \(\(:\)\(\ \)\), "\<\"Possible spelling error: new symbol name \\\"\\!\\(Len\\)\\\" is similar to existing symbol \\\"\\!\\(len\\)\\\". \\!\\(\\*ButtonBox[\\\"More\[Ellipsis]\\\", ButtonStyle->\\\"RefGuideLinkText\\\", ButtonFrame->None, ButtonData:>\\\"General::spell1\\\"]\\)\"\>"}]\)
*)
ShowSpinEigenfunctionCompatibility[] := Norsym
(*
NN is the number of electrons
*)
(*
\[Theta] is a list of all possible eigenfunction of \!\(TraditionalForm\`S\_z\) with a given number of spin up and down components.  The number of these \
functions is given by PermutationNumber.
*)
(*
MM is a matrix whose eigenvectors give eigenfunctions of \!\(TraditionalForm\`S\^2\).  It is calculated using the Dirac Identity which expresses the \!\(TraditionalForm\`S\^2\)operator in terms of permutations.
*)
(*
VV gives a list of orthonormal Spin Eigenfunctions with the total spin \
quantum number S.
*)
(*
PP gives a list of all possible permutations on NN electrons with ordered \
with alternating parity.  There are NN! of these. 
*)
(*
The above functions transform the spin eigenfunction basis to one in which \
each function is either symmetric or antisymmetric when the electrons 2i and \
2i-1 are transposed i = 1... \!\(TraditionalForm\`Floor[NN\/2]\).  This is possible as these permutations commute.  
*)
(*
Norsym is a array that tells us which spin eigenfunctions vanish when the \
spatial function has one or more doubly occupied orbitals.  The \!\(TraditionalForm\`n\^\(\(th\)\(\ \)\)\)table corresponds to having n pairs of doubly occupied orbitals.  The \!\(TraditionalForm\`i\^th\)entry corresponds to the \!\(TraditionalForm\`i\^th\)spin eigenfunction and if it is zero then that eigenfunction vanishes.
*)
\!\(\(\(\[IndentingNewLine]\)\(SpinVectorConvert[v_List]\  := Module[{}, \ \[Theta]\[Theta]\  = \((\ \[Theta]\  /. \ {1 -> \ \[Alpha], \ 0 -> \ \[Beta]})\); \[IndentingNewLine]\[Chi] = \ {}; \[IndentingNewLine]For[i = 1, \ i <= \ PermutationNumber, \ \ \[Chi]\  = \ Append[\[Chi], \ ReplacePart[\[Theta]\[Theta][\([i]\)], AngleBracket, 0]]; \ \(i++\)]; \[IndentingNewLine]\ \ \[Sum]\+\(i = 1\)\%PermutationNumber v[\([i]\)]\ \[Chi][\([i]\)]]\)\)\)
\!\(\(\(\[Psi]\_i_[VV_] := \ \[Sum]\+\(jj = 1\)\%\(Length[VV]\)\ \((SpinVectorConvert[VV[\([jj]\)]]\ \@\(1\/\(NN!\)\)\ \(\[Sum]\+\(p = 1\)\%\(NN\ !\)\((\(-1\))\)\^\(p - 1\)\ \(UUU[p, VV]\)[\([jj, i]\)]\ ReplacePart[PP[\([p]\)], \ BracketingBar, 0]\))\)\)\(\[IndentingNewLine]\)
  \)\)
(*
SpinVectorConvert converts the representation of a Spin Eigenfunction in \
terms of a vector into a more explicit form in terms of a linear combination \
of basic eigenfunctions of \!\(TraditionalForm\`S\_z\).
\!\(\*FormBox[
  StyleBox[\(\[Psi]\_i\),
    FontWeight->"Bold"], TraditionalForm]\) gives the antisymmetrized wavefunction corresponding to the \!\(TraditionalForm\`i\^th\)spin eigenfunction.  These are not normalized for double occupied orbitals.
*)
(*
Norsym is a array that tells us which spin eigenfunctions vanish when the \
spatial function has one or more doubly occupied orbitals.  The \!\(TraditionalForm\`n\^\(\(th\)\(\ \)\)\)table corresponds to having n pairs of doubly occupied orbitals.  The \!\(TraditionalForm\`i\^th\)entry corresponds to the \!\(TraditionalForm\`i\^th\)spin eigenfunction and if it is zero then that eigenfunction vanishes.
*)
ShowSpinBasis[] := 
  Table[SpinVectorConvert[Vtrans[[n]]], {n, 1, Length[Vtrans]}] // TableForm
(*
We can see the calculated Spin Eigenfunctions above.  The \[Alpha] and \
\[Beta] in represent spin up and spin down electron.  We can also notice that \
all are either symmetric or antisymmetric under the exchange of the first and \
second electron and also the third and fourth electron, etc.
*)
RepresentationMatrix[perm_List] := UUU[Position[PP, perm][[1, 1]] , Vtrans]
(*
3.8 File Saving
*)
(*
This section saves intermediate data so that future calculations can be \
accelerated.
*)
SaveIntegralandHamiltonianData := Module[{},
       Definition[Idef] >> Ideffiletemp;
      If[FileByteCount["Ideffile"] < FileByteCount["Ideffiletemp"], 
        DeleteFile["Ideffile"]; RenameFile["Ideffiletemp", "Ideffile"]];
      
      alfexact >> alfadataexacttemp;
      If[FileByteCount["alfadataexact"] < FileByteCount["alfadataexacttemp"], 
        DeleteFile["alfadataexact"]; 
        RenameFile["alfadataexacttemp", "alfadataexact"]];
      ];
(*
End
*)
End[]
(*
SACI`Private`
*)
EndPackage[]
