(* Content-type: application/vnd.wolfram.mathematica *)

(*** Wolfram Notebook File ***)
(* http://www.wolfram.com/nb *)

(* CreatedBy='Mathematica 9.0' *)

(*CacheID: 234*)
(* Internal cache information:
NotebookFileLineBreakTest
NotebookFileLineBreakTest
NotebookDataPosition[       157,          7]
NotebookDataLength[     15663,        401]
NotebookOptionsPosition[     14531,        361]
NotebookOutlinePosition[     14975,        378]
CellTagsIndexPosition[     14932,        375]
WindowFrame->Normal*)

(* Beginning of Notebook Content *)
Notebook[{

Cell[CellGroupData[{
Cell["\<\
Calculate energies of three lowest states of the He atom using 3 basis states \
(energies are in atomic unit)\
\>", "Title",
 CellChangeTimes->{{3.5720155839971457`*^9, 3.572015596770288*^9}, {
  3.594837777572315*^9, 3.594837799695654*^9}, {3.7672534115080605`*^9, 
  3.7672534437065163`*^9}, {3.76725374977905*^9, 3.7672537666582794`*^9}, {
  3.767255654386371*^9, 3.7672556665387926`*^9}}],

Cell["\<\
Note: we show the codes for the calculations using 3 basis states here so \
that the result can be directly compared to our analytic calculation in our \
research article. Line 2 of the code can then be changed to various number of \
basis states. In our article, we vary the number of basis states from n = 3 \
to n = 25.

Define radial wave function, i.e. the equation (16) in the research article.\
\>", "Text",
 CellChangeTimes->{{3.57201560320553*^9, 3.572015612216217*^9}, {
  3.572016598152993*^9, 3.572016600484062*^9}, {3.57201718859979*^9, 
  3.572017190328268*^9}, {3.7672534574501405`*^9, 3.7672534908029985`*^9}, {
  3.7672541249909043`*^9, 3.76725412815771*^9}, {3.767254923805897*^9, 
  3.7672550502129173`*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{
   SubscriptBox["P", "n_"], "[", "r_", "]"}], ":=", 
  RowBox[{
   SqrtBox[
    FractionBox[
     RowBox[{"2", 
      RowBox[{
       RowBox[{"(", 
        RowBox[{"n", "-", "1"}], ")"}], "!"}], " "}], 
     RowBox[{
      SuperscriptBox["n", "2"], 
      RowBox[{
       RowBox[{"(", "n", ")"}], "!"}]}]]], 
   RowBox[{"(", 
    FractionBox[
     RowBox[{"4", " ", "r"}], "n"], ")"}], 
   SuperscriptBox["\[ExponentialE]", 
    RowBox[{"-", 
     FractionBox[
      RowBox[{"2", "r"}], "n"]}]], 
   RowBox[{"LaguerreL", "[", 
    RowBox[{
     RowBox[{"n", "-", "1"}], ",", "1", ",", 
     FractionBox[
      RowBox[{"4", "r"}], "n"]}], "]"}]}]}]], "Input",
 CellChangeTimes->{{3.571682606534275*^9, 3.571682808440976*^9}, {
   3.571683051432061*^9, 3.571683053639886*^9}, {3.5716831963545227`*^9, 
   3.571683206297831*^9}, {3.571683264033547*^9, 3.571683264439849*^9}, {
   3.571684711395982*^9, 3.571684728792845*^9}, {3.571684813835341*^9, 
   3.57168481672513*^9}, 3.571771419917448*^9, 3.572014882532148*^9, {
   3.6345555711955643`*^9, 3.634555574869687*^9}, {3.767253493876204*^9, 
   3.767253548257899*^9}}],

Cell["\<\
Generate the list of states. Here, we use 3 basis states as shown in equation \
(19) in the article. This can be changed to various number of basis states. \
\>", "Text",
 CellChangeTimes->{{3.586111549982008*^9, 3.5861115567561626`*^9}, {
  3.595556196859631*^9, 3.59555620534161*^9}, {3.5956775574322777`*^9, 
  3.595677560255501*^9}, {3.595679557833642*^9, 3.595679702789641*^9}, {
  3.767254113368884*^9, 3.7672541143672857`*^9}, {3.7672541455361404`*^9, 
  3.767254186704612*^9}, {3.767254246359116*^9, 3.7672542720523605`*^9}, {
  3.7672547466519885`*^9, 3.76725474744759*^9}, {3.7672557766777887`*^9, 
  3.7672557790801926`*^9}}],

Cell[CellGroupData[{

Cell[BoxData[
 RowBox[{"st", "=", 
  RowBox[{"DeleteDuplicates", "[", 
   RowBox[{
    RowBox[{"Table", "[", 
     RowBox[{
      RowBox[{"{", 
       RowBox[{
        RowBox[{"{", 
         RowBox[{"1", ",", "i"}], "}"}], ",", 
        RowBox[{"{", 
         RowBox[{"i", ",", "1"}], "}"}]}], "}"}], ",", 
      RowBox[{"{", 
       RowBox[{"i", ",", "2"}], "}"}]}], "]"}], "~", "Flatten", "~", "1"}], 
   "]"}]}]], "Input",
 CellChangeTimes->{{3.5716838814593163`*^9, 3.571684015453298*^9}, {
   3.571684046035569*^9, 3.571684058832741*^9}, {3.5716841230899754`*^9, 
   3.571684376504847*^9}, {3.571684441190647*^9, 3.5716844451504107`*^9}, 
   3.571684554709581*^9, {3.571684588210924*^9, 3.571684609391206*^9}, 
   3.5716868451044807`*^9, {3.5716875302649717`*^9, 3.5716875390818777`*^9}, {
   3.57168769091436*^9, 3.571687721674882*^9}, {3.5716877691268797`*^9, 
   3.571687774936057*^9}, {3.571687900771812*^9, 3.57168790116206*^9}, {
   3.571688001546517*^9, 3.571688003110647*^9}, {3.5716880368969793`*^9, 
   3.5716880469433537`*^9}, {3.571751745126295*^9, 3.571751758864758*^9}, {
   3.571760010237328*^9, 3.5717600491333427`*^9}, 3.5717600961733847`*^9, {
   3.5717716191262703`*^9, 3.571771631493891*^9}, {3.57177166188062*^9, 
   3.571771698528517*^9}, {3.5717717566063643`*^9, 3.571771757178018*^9}, {
   3.572014895174371*^9, 3.5720149282425547`*^9}, {3.572015246762957*^9, 
   3.572015300492633*^9}, {3.572015378766746*^9, 3.572015378906934*^9}, {
   3.5720154193526077`*^9, 3.5720154426042013`*^9}, {3.57201550851149*^9, 
   3.572015508780752*^9}, 3.5720933344386787`*^9, {3.586092721852374*^9, 
   3.5860927341225147`*^9}, {3.586092934924246*^9, 3.586092946810019*^9}, {
   3.587201332531014*^9, 3.587201350580724*^9}, {3.594837677485523*^9, 
   3.594837681163014*^9}, {3.594837926222932*^9, 3.594837942102412*^9}, {
   3.59483797373253*^9, 3.594837995830896*^9}, {3.594838083448278*^9, 
   3.594838126039194*^9}, 3.594839454888926*^9, 3.594839534263056*^9, {
   3.594839652713702*^9, 3.594839746809248*^9}, {3.594840150269005*^9, 
   3.594840154000621*^9}, {3.5955560868169947`*^9, 3.5955560878811493`*^9}, {
   3.6350991339818363`*^9, 3.635099134282207*^9}, {3.7672535875387673`*^9, 
   3.7672536021715927`*^9}, 3.767253884953286*^9, 3.767253920677348*^9, 
   3.7672540736512146`*^9, 3.7672542232710757`*^9, {3.7672550805861707`*^9, 
   3.76725508613978*^9}, {3.7672551259666495`*^9, 3.7672551416602774`*^9}, 
   3.7672553396714225`*^9}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"{", 
    RowBox[{"1", ",", "1"}], "}"}], ",", 
   RowBox[{"{", 
    RowBox[{"1", ",", "2"}], "}"}], ",", 
   RowBox[{"{", 
    RowBox[{"2", ",", "1"}], "}"}]}], "}"}]], "Output",
 CellChangeTimes->{3.767254874884212*^9, 3.767255243466055*^9, 
  3.7672553825246973`*^9, 3.7672556051994853`*^9}]
}, Open  ]],

Cell["\<\
Calculate unperturbed energy of the He+ ion, i.e. equation (11) in the text\
\>", "Text",
 CellChangeTimes->{{3.594840023888061*^9, 3.594840025728023*^9}, {
  3.594840180095477*^9, 3.594840182966941*^9}, {3.595678003541584*^9, 
  3.595678003731388*^9}, {3.5956791527994213`*^9, 3.595679166264566*^9}, {
  3.595679198011915*^9, 3.5956793104897747`*^9}, {3.595679392874557*^9, 
  3.59567949631775*^9}, {3.5956797150183687`*^9, 3.595679716728313*^9}, {
  3.767254302488014*^9, 3.7672543518933*^9}, {3.7672547369019713`*^9, 
  3.7672547376195726`*^9}, {3.7672548194417152`*^9, 3.7672548379901476`*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{"E0", "[", 
   RowBox[{"{", 
    RowBox[{"n1_", ",", "n2_"}], "}"}], "]"}], ":=", 
  RowBox[{
   FractionBox[
    RowBox[{"-", "2"}], 
    SuperscriptBox["n1", "2"]], "+", 
   FractionBox[
    RowBox[{"-", "2"}], 
    SuperscriptBox["n2", "2"]]}]}]], "Input",
 CellChangeTimes->{{3.5948382270417356`*^9, 3.59483825261222*^9}, {
   3.594839016978347*^9, 3.59483904302036*^9}, {3.594841257090726*^9, 
   3.59484128564902*^9}, {3.5948422643479*^9, 3.5948423431129436`*^9}, 
   3.7672542381223016`*^9}],

Cell["\<\
Evaluate the matrix elements of the potential, i.e. the second term of the \
equation (14) in the text. \
\>", "Text",
 CellChangeTimes->{{3.594840052549675*^9, 3.594840071054782*^9}, 
   3.594841301973468*^9, {3.595679743986805*^9, 3.595679749382815*^9}, {
   3.767254378709747*^9, 3.7672544261338296`*^9}, {3.7672547285091567`*^9, 
   3.7672547290083575`*^9}}],

Cell[BoxData[
 RowBox[{
  RowBox[{
   RowBox[{"V", "[", 
    RowBox[{
     RowBox[{"{", 
      RowBox[{"n1_", ",", "n2_"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"n3_", ",", "n4_"}], "}"}]}], "]"}], ":=", 
   "\[IndentingNewLine]", 
   RowBox[{"NIntegrate", "[", 
    RowBox[{
     RowBox[{
      RowBox[{
       SubscriptBox["P", "n1"], "[", "r1", "]"}], 
      RowBox[{
       SubscriptBox["P", "n3"], "[", "r1", "]"}], 
      RowBox[{"Min", "[", 
       RowBox[{
        FractionBox["1", "r2"], ",", 
        FractionBox["1", "r1"]}], "]"}], 
      RowBox[{
       SubscriptBox["P", "n2"], "[", "r2", "]"}], 
      RowBox[{
       SubscriptBox["P", "n4"], "[", "r2", "]"}]}], ",", 
     RowBox[{"{", 
      RowBox[{"r2", ",", "0", ",", "\[Infinity]"}], "}"}], ",", 
     RowBox[{"{", 
      RowBox[{"r1", ",", "0", ",", "\[Infinity]"}], "}"}]}], "]"}]}], 
  ";"}]], "Input",
 CellChangeTimes->{{3.594838876557716*^9, 3.5948389943276453`*^9}, {
   3.594839067133453*^9, 3.594839092550424*^9}, {3.594839391254321*^9, 
   3.594839391982773*^9}, 3.594839571778838*^9, {3.5948397800892277`*^9, 
   3.59483978216215*^9}, {3.594841308965118*^9, 3.59484133370759*^9}, {
   3.7672538698056593`*^9, 3.7672538919888983`*^9}, 3.7672542330366926`*^9, 
   3.767254347353692*^9, {3.7672551044698124`*^9, 3.767255106700616*^9}, 
   3.767255159335108*^9}],

Cell["\<\
Calculate all the elements of the Hamiltonian matrix, i.e. equation (14) in \
the text\
\>", "Text",
 CellChangeTimes->{{3.595677566094204*^9, 3.5956775695255203`*^9}, {
   3.595680121009981*^9, 3.595680204178039*^9}, {3.7672544457586637`*^9, 
   3.767254519297192*^9}, {3.7672546737998614`*^9, 3.7672546743770623`*^9}, {
   3.767254786666058*^9, 3.767254812920904*^9}, {3.767255411337947*^9, 
   3.767255415300354*^9}, 3.767255524687745*^9}],

Cell[BoxData[
 RowBox[{
  RowBox[{"H", "=", 
   RowBox[{
    RowBox[{"DiagonalMatrix", "[", 
     RowBox[{"E0", "/@", "st"}], "]"}], "+", 
    RowBox[{"Table", "[", 
     RowBox[{
      RowBox[{"V", "[", 
       RowBox[{
        RowBox[{"st", "\[LeftDoubleBracket]", "i", "\[RightDoubleBracket]"}], 
        ",", 
        RowBox[{
        "st", "\[LeftDoubleBracket]", "j", "\[RightDoubleBracket]"}]}], "]"}],
       ",", 
      RowBox[{"{", 
       RowBox[{"i", ",", 
        RowBox[{"Length", "[", "st", "]"}]}], "}"}], ",", 
      RowBox[{"{", 
       RowBox[{"j", ",", 
        RowBox[{"Length", "[", "st", "]"}]}], "}"}]}], "]"}]}]}], 
  ";"}]], "Input",
 CellChangeTimes->{{3.5948413620696297`*^9, 3.594841365920395*^9}, {
   3.594841698295514*^9, 3.5948417025888023`*^9}, {3.594841744982246*^9, 
   3.5948417725429277`*^9}, {3.595553162705166*^9, 3.595553163464281*^9}, 
   3.7672538660772533`*^9, 3.767253922876952*^9, 3.7672543453412886`*^9, 
   3.7672544426386585`*^9, 3.7672551609419107`*^9, 3.7672551931403666`*^9}],

Cell["\<\
Show the Hamiltonian matrix which can be compared to analytic calculations in \
equation (29) in the text\
\>", "Text",
 CellChangeTimes->{{3.595677574306884*^9, 3.5956775769055557`*^9}, {
  3.595677649536332*^9, 3.595677786223404*^9}, {3.595680050877356*^9, 
  3.5956800710590773`*^9}, {3.767254575036089*^9, 3.7672546275145807`*^9}, {
  3.7672547575096073`*^9, 3.7672547835928526`*^9}, {3.7672554063147383`*^9, 
  3.767255408935543*^9}, {3.767255531021356*^9, 3.767255558040603*^9}}],

Cell[CellGroupData[{

Cell[BoxData[
 RowBox[{"H", "//", "MatrixForm"}]], "Input",
 CellChangeTimes->{{3.594841450620315*^9, 3.594841455525065*^9}, {
   3.5948416120951*^9, 3.594841626931789*^9}, {3.594841774678866*^9, 
   3.594841780785747*^9}, {3.59484193819151*^9, 3.594842086969365*^9}, {
   3.594842121175425*^9, 3.594842121506407*^9}, {3.595679911313449*^9, 
   3.595679970700964*^9}, 3.595680046975935*^9, {3.595680097208321*^9, 
   3.595680099108556*^9}, 3.639245340137539*^9, 3.767253860788844*^9, 
   3.7672545643344707`*^9, {3.7672555905042596`*^9, 3.7672555957146688`*^9}}],

Cell[BoxData[
 TagBox[
  RowBox[{"(", "\[NoBreak]", GridBox[{
     {
      RowBox[{"-", "2.749999566994383`"}], "0.1787100634316854`", 
      "0.17871006666783404`"},
     {"0.1787100634316854`", 
      RowBox[{"-", "2.080246912964298`"}], "0.043895746665815756`"},
     {"0.17871006666783404`", "0.043895746665815756`", 
      RowBox[{"-", "2.0802469126922434`"}]}
    },
    GridBoxAlignment->{
     "Columns" -> {{Center}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, 
      "RowsIndexed" -> {}},
    GridBoxSpacings->{"Columns" -> {
        Offset[0.27999999999999997`], {
         Offset[0.7]}, 
        Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
        Offset[0.2], {
         Offset[0.4]}, 
        Offset[0.2]}, "RowsIndexed" -> {}}], "\[NoBreak]", ")"}],
  Function[BoxForm`e$, 
   MatrixForm[BoxForm`e$]]]], "Output",
 CellChangeTimes->{
  3.594841629513722*^9, 3.5948417174655*^9, {3.594841782541287*^9, 
   3.594841789120206*^9}, {3.594842101317523*^9, 3.594842124526231*^9}, 
   3.5948423882656107`*^9, 3.595553117145767*^9, 3.595554298276093*^9, 
   3.595679788564268*^9, 3.6350991447185707`*^9, 3.639245347423904*^9, 
   3.767254878425418*^9, 3.767255248239663*^9, 3.767255386299904*^9, 
   3.767255608553491*^9}]
}, Open  ]],

Cell[TextData[{
 "Calculate and show the eigenvalues for 1s1s, ",
 Cell[BoxData[
  FormBox[
   RowBox[{Cell["1s2s"], " ", "triplet"}], TraditionalForm]]],
 ", ",
 Cell[BoxData[
  FormBox[
   RowBox[{"and", "  ", Cell["1s2s"], " ", "singlet"}], TraditionalForm]]],
 "; the results are then compared to the analytic caculations shown in \
equation (30) in the text"
}], "Text",
 CellChangeTimes->{{3.595677832584819*^9, 3.595677836567211*^9}, 
   3.7672546676534505`*^9, 3.767255484408475*^9}],

Cell[CellGroupData[{

Cell[BoxData[
 RowBox[{
  RowBox[{"eval", "=", 
   RowBox[{
    RowBox[{"Eigenvalues", "[", "H", "]"}], "\[LeftDoubleBracket]", 
    RowBox[{";;", "3"}], "\[RightDoubleBracket]"}]}], ";", "eval"}]], "Input",\

 CellChangeTimes->{{3.594841450620315*^9, 3.594841455525065*^9}, {
   3.5948416120951*^9, 3.594841626931789*^9}, {3.594841774678866*^9, 
   3.594841780785747*^9}, {3.59484193819151*^9, 3.594842086969365*^9}, {
   3.594842121175425*^9, 3.594842121506407*^9}, {3.595679911313449*^9, 
   3.595679970700964*^9}, 3.595680046975935*^9, {3.595680097208321*^9, 
   3.595680099108556*^9}, 3.639245340137539*^9, 3.767253860788844*^9, 
   3.7672545643344707`*^9}],

Cell[BoxData[
 RowBox[{"{", 
  RowBox[{
   RowBox[{"-", "2.8304373986629234`"}], ",", 
   RowBox[{"-", "2.124142659494087`"}], ",", 
   RowBox[{"-", "1.9559133344939137`"}]}], "}"}]], "Output",
 CellChangeTimes->{3.7672556086002913`*^9}]
}, Open  ]]
}, Open  ]]
},
WindowSize->{1213, 710},
WindowMargins->{{Automatic, 1}, {Automatic, 51}},
ShowSelection->True,
FrontEndVersion->"10.3 for Microsoft Windows (32-bit) (December 10, 2015)",
StyleDefinitions->FrontEnd`FileName[{"Report"}, "StandardReport.nb", 
  CharacterEncoding -> "WindowsANSI"]
]
(* End of Notebook Content *)

(* Internal cache information *)
(*CellTagsOutline
CellTagsIndex->{}
*)
(*CellTagsIndex
CellTagsIndex->{}
*)
(*NotebookFileOutline
Notebook[{
Cell[CellGroupData[{
Cell[579, 22, 402, 7, 164, "Title"],
Cell[984, 31, 737, 13, 87, "Text"],
Cell[1724, 46, 1152, 33, 66, "Input"],
Cell[2879, 81, 646, 10, 30, "Text"],
Cell[CellGroupData[{
Cell[3550, 95, 2457, 42, 40, "Input"],
Cell[6010, 139, 345, 10, 40, "Output"]
}, Open  ]],
Cell[6370, 152, 608, 9, 30, "Text"],
Cell[6981, 163, 532, 15, 56, "Input"],
Cell[7516, 180, 372, 7, 30, "Text"],
Cell[7891, 189, 1349, 36, 76, "Input"],
Cell[9243, 227, 452, 8, 30, "Text"],
Cell[9698, 237, 1027, 26, 40, "Input"],
Cell[10728, 265, 495, 8, 30, "Text"],
Cell[CellGroupData[{
Cell[11248, 277, 562, 8, 40, "Input"],
Cell[11813, 287, 1257, 29, 66, "Output"]
}, Open  ]],
Cell[13085, 319, 491, 13, 30, "Text"],
Cell[CellGroupData[{
Cell[13601, 336, 662, 13, 40, "Input"],
Cell[14266, 351, 237, 6, 40, "Output"]
}, Open  ]]
}, Open  ]]
}
]
*)

(* End of internal cache information *)