Hydrogen Atom Quantum Mechanics Done With Optimization

Posted 1 month ago
250 Views
|
0 Replies
|
3 Total Likes
|
 In[5]:= sphcmpt[pefunc, nInteger, d_] := Block[{i, vars, psi, vwb, ener, enres, psires, varres, cons1, cons2, sln, wavefunctions}, (* define variables *) vars = Table[psi[r], {r, 0, 1, d}];  (* add bounds to variables *) vwb = {#, -1000, 1000} & /@ vars; (* calculate expectation values of energy for a radial function \ using finite differences with spacing d for derivatives *) ener = d Sum[-psi[ r] (2/r ((psi[r + d] - psi[r - d])/(2 d)) + (( psi[r + d] - 2*psi[r] + psi[r - d])/d^2)) 4 \[Pi] r^2 + psi[r] pefunc psi[r] 4 \[Pi] r^2, {r, d, 1 - d, d}]; (* allocate variables for results *) enres = Table[0, n]; wavefunctions = Table[0, n]; (* boundary condition and normalization condition *) cons1 = {psi[1] == 0, 4 \[Pi] Sum[r^2 psi[r]^2, {r, d, 1, d}] == d^-1}; (* calculate lowest n energy levels *) Do[ (* orthogonality constraints *) cons2 = Table[Sum[ psi[r]*(psi[r] /. wavefunctions[[j]])*r^2, {r, 0, 1, d}] == 0, {j, i - 1}]; (* use ParametricIPOPTMinimize to solve problem, using 10 random searches, seed 0 *) sln = iMin[ener, Join[cons1, cons2], vwb, 10, 0]; (* save results *) enres[[i]] = sln[[1]]; wavefunctions[[i]] = sln[[2]], {i, n} ]; {enres, wavefunctions} ] In[6]:= (* perform calculation *) hatom = sphcmpt[(-100)/r, 5, 10^-3]; During evaluation of In[6]:= {{Solve_Succeeded,10}} During evaluation of In[6]:= {{Solve_Succeeded,10}} During evaluation of In[6]:= {{Solve_Succeeded,10}} During evaluation of In[6]:= {{Solve_Succeeded,10}} During evaluation of In[6]:= {{Solve_Succeeded,10}} In[7]:= (* look at energy levels *) hatom[[1]] Out[7]= {-2498.44, -624.902, -277.758, -156.015, -89.3199} In[8]:= (* check that energy levels go as 1/n^2 *) Table[hatom[[1, i]]*i^2, {i, 5}] Out[8]= {-2498.44, -2499.61, -2499.83, -2496.25, -2233.} In[9]:= (* Plot wavefunctions *) Table[ListPlot[Table[{r, psi[r]}, {r, 0, 1, 1/100}] /. hatom[[2, n]], Joined -> True, PlotRange -> All], {n, 5}]  (* Verify that lowest state is exponential in r *) ListPlot[Table[{r, Log @ psi[r]}, {r, 0, 1, 1/100}] /. hatom[[2, 1]]]