The atomic orbitals are solutions to the equation,

\begin{equation} - \frac{\hbar^2}{2m}\nabla^2\phi(\vec{r}) - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}|}\phi(\vec{r}) = E\phi(\vec{r}). \end{equation}

In spherical coordinates the Laplacian operator is,

\[\begin{equation} \nabla ^2 \phi = \frac{1}{r^2} \frac{\partial}{\partial r}\left(r^2 \frac{\partial \phi}{\partial r}\right)+\frac{1}{r^2 \sin\theta} \frac{\partial}{\partial \theta}\left(\sin \theta~\frac{\partial \phi}{\partial \theta}\right)+\frac{1}{r^2 \sin^2\theta} \frac{\partial^2 \phi}{\partial \varphi^2} \end{equation}\]

Often the atomic orbitials are used to evaluate the matrix elements of the form, $\langle \phi_m |H|\phi_n\rangle$. To perform this calculation it is necessary to calculate $\nabla^2 \phi_n$. Since numerical differentiation can be imprecise, it is useful to calculate the Laplacian analytically. To facilitate such calculations, a table of atomic orbitals and the Laplacian of the orbitals is given below.

\[\begin{equation} \begin{array}{c,c} \phi_{1s}=\sqrt{\frac{Z^3}{\pi a^3_0}}\exp\left(-\frac{Zr}{a_0}\right)& \nabla ^2\phi_{1s}=\sqrt{\frac{Z^5}{\pi a^5_0}} \left(\frac{Z}{a_0}-\frac{2}{r}\right)\exp\left(-\frac{Zr}{a_0}\right)\\ \phi_{2s}=\frac{1}{4} \sqrt{\frac{Z^3}{2\pi a_0^3}}\left(2-\frac{Zr}{a_0}\right)\exp\left(-\frac{Zr}{2a_0}\right)& \nabla ^2\phi_{2s}=\frac{1}{4}\sqrt{\frac{Z^5}{2\pi a^5_0}} \left(-\frac{Z^2r}{4a_0^2}+\frac{5Z}{2a_0}-\frac{4}{r}\right)\exp\left(-\frac{Zr}{2a_0}\right)\\ \phi_{2p_x}=\frac{1}{4} \sqrt{\frac{Z^5}{2\pi a_0^5}}r \exp\left(-\frac{Zr}{2a_0}\right) \sin\theta \cos\varphi & \nabla^2\phi_{2p_x}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\sin\theta\cos\varphi \exp\left(-\frac{Zr}{2a_0}\right) \\ \phi_{2p_y}=\frac{1}{4} \sqrt{\frac{Z^5}{2\pi a_0^5}}r \exp\left(-\frac{Zr}{2a_0}\right) \sin\theta \sin\varphi & \nabla^2\phi_{2p_y}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\sin\theta\sin\varphi \exp\left(-\frac{Zr}{2a_0}\right)\\ \phi_{2p_z}=\frac{1}{4} \sqrt{\frac{Z^5}{2\pi a_0^5}}r \exp\left(-\frac{Zr}{2a_0}\right) \cos\theta& \nabla^2\phi_{2p_z}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\cos\theta \exp\left(-\frac{Zr}{2a_0}\right)\\ \phi_{3s} = \frac{1}{81}\sqrt{\frac{Z^3}{3\pi a_0^3}} \left(27-\frac{18Zr}{a_0}+\frac{2Z^2r^2}{a_0^2}\right) \exp\left(\frac{-Zr}{3a_0}\right) & \nabla^2\phi_{3s} = \frac{1}{81}\sqrt{\frac{Z^3}{3\pi a_0^3}} \left(-\frac{54 Z}{a_0 r}+\frac{39Z^2}{a_0^2}-\frac{6Z^3r}{a_0^3}+\frac{2Z^4r^2}{9a_0^4}\right) \exp\left(\frac{-Zr}{3a_0}\right) \\ \phi_{3p_x} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(6-\frac{Zr}{a_0}\right)\frac{Zr}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \sin\theta\cos\varphi & \nabla^2\phi_{3p_x} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(-\frac{12Z}{a_0r}+\frac{8Z^2}{3a_0^2}-\frac{Z^3r}{9a_0^3}\right) \frac{Zr}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \sin\theta\cos\varphi \\ \phi_{3p_y} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(6-\frac{Zr}{a_0}\right)\frac{Zr}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \sin\theta\sin\varphi & \nabla^2\phi_{3p_y} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(-\frac{12Z^2}{a_0^2}+\frac{8Z^3r}{3a_0^3}-\frac{Z^4r^2}{9a_0^4}\right) \exp\left(\frac{-Zr}{3a_0}\right) \sin\theta\sin\varphi \\ \phi_{3p_z} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(6-\frac{Zr}{a_0}\right)\frac{Zr}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \cos\theta & \nabla^2\phi_{3p_z} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(-\frac{12Z^2}{a_0^2}+\frac{8Z^3r}{3a_0^3}-\frac{Z^4r^2}{9a_0^4}\right) \exp\left(\frac{-Zr}{3a_0}\right) \cos\theta \\ \phi_{3d_{z^2}} = \frac{1}{81}\sqrt{\frac{Z^7}{6\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\left(3\cos^2\left(\theta\right)-1\right) & \nabla^2\phi_{3d_{z^2}} = \frac{1}{243}\sqrt{\frac{Z^9}{6\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\left(3\text{cos}^2\left(\theta\right)-1\right)\\ \phi_{3d_{xz}} = \frac{1}{81}\sqrt{\frac{2Z^7}{\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\sin\left(\theta\right)\cos\left(\theta\right)\cos\left(\phi\right) & \nabla^2\phi_{3d_{xz}} = \frac{1}{243}\sqrt{\frac{2Z^9}{\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\sin\left(\theta\right)\cos\left(\theta\right)\cos\left(\phi\right)\\ \phi_{3d_{yz}} = \frac{1}{81}\sqrt{\frac{2Z^7}{\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\sin\left(\theta\right)\cos\left(\theta\right)\sin\left(\phi\right) & \nabla^2\phi_{3d_{yz}} = \frac{1}{243}\sqrt{\frac{2Z^9}{\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\sin\left(\theta\right)\cos\left(\theta\right)\sin\left(\phi\right)\\ \phi_{3d_{xy}} = \frac{1}{81}\sqrt{\frac{2Z^7}{\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\sin^2\left(\theta\right)\sin\left(\phi\right)\cos\left(\phi\right) & \nabla^2\phi_{3d_{xy}} = \frac{1}{243}\sqrt{\frac{2Z^9}{\pi a_0^9}}~r\left(\frac{Zr}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\sin^2\left(\theta\right)\sin\left(\phi\right)\cos\left(\phi\right)\\ \phi_{3d_{x^2-y^2}} = \frac{1}{81}\sqrt{\frac{Z^7}{2\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\sin^2\left(\theta\right)\text{cos}\left(2\phi\right) & \nabla^2\phi_{3d_{x^2-y^2}} = \frac{1}{243}\sqrt{\frac{Z^9}{2\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\sin^2\left(\theta\right)\cos\left(2\phi\right)\\ \end{array} \end{equation}\]

A project has been established to collect some code for calculating with atomic orbitals. There are programs to check that the atomic orbitals are defined correctly for the 1s and 2p_z atomic orbitals. In the code block below, the atomic orbitals $\phi (r,\theta ,\varphi )$ are defined and then 10 random positions are chosen for $r$, $\theta$, and $\varphi$. The energy $E=H\phi/\phi$ is calculated at these positions. If the wave functions are properly programmed, the energy should be the same at all of these positions.

If the JavaScript code on the left is modified, pressing the Execute button displays the results of executing the modified code.

Only JavaScript will run in a browser. Copy the code on the left and run it in an appropriate environment. <script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = 3.1415926535897932; var hbar = 1.05457266E-34; var ec = 1.60217733E-19; var eps0 = 8.854187817E-12; var me = 9.1093897E-31; var k1 = 0.5*hbar*hbar/me; var k2 = ec*ec/(4*pi*eps0); //Coulomb constant function ao_1s(r,theta,phi,Z) { //the 1s atomic orbital return Math.sqrt(Math.pow(Z,3)/(pi*Math.pow(a0,3)))*Math.exp(-Z*r/a0); } function L_ao_1s(r,theta,phi,Z) { //Lagrangian of the 1s atomic orbital return Math.sqrt(Math.pow(Z,5)/(pi*Math.pow(a0,5)))*(Z/a0 -2/r)*Math.exp(-Z*r/a0); } function ao_2s(r,theta,phi,Z) { //the 2s atomic orbital return 0.25*Math.sqrt(Math.pow(Z,3)/(2*pi*Math.pow(a0,3)))*(2-Z*r/a0)*Math.exp(-Z*r/(2*a0)); } function L_ao_2s(r,theta,phi,Z) { //Lagrangian of the 2s atomic orbital return 0.25*Math.sqrt(Math.pow(Z,5)/(2*pi*Math.pow(a0,5)))*(-Z*Z*r/(4*a0*a0)+5*Z/(2*a0)-4/r)*Math.exp(-Z*r/(2*a0)); } function ao_2px(r,theta,phi,Z) { //the 2px atomic orbital return 0.25*Math.sqrt(Math.pow(Z,5)/(2*pi*Math.pow(a0,5)))*r*Math.sin(theta)*Math.cos(phi)*Math.exp(-Z*r/(2*a0)); } function L_ao_2px(r,theta,phi,Z) { //Lagrangian of the 2px atomic orbital return Math.sqrt(Math.pow(Z,7)/(32*pi*Math.pow(a0,7)))*(Z*r/(4*a0)-2)*Math.sin(theta)*Math.cos(phi)*Math.exp(-Z*r/(2*a0)); } function ao_2py(r,theta,phi,Z) { //the 2py atomic orbital return 0.25*Math.sqrt(Math.pow(Z,5)/(2*pi*Math.pow(a0,5)))*r*Math.sin(theta)*Math.sin(phi)*Math.exp(-Z*r/(2*a0)); } function L_ao_2py(r,theta,phi,Z) { //Lagrangian of the 2py atomic orbital return Math.sqrt(Math.pow(Z,7)/(32*pi*Math.pow(a0,7)))*(Z*r/(4*a0)-2)*Math.sin(theta)*Math.sin(phi)*Math.exp(-Z*r/(2*a0)); } function ao_2pz(r,theta,phi,Z) { //the 2pz atomic orbital return 0.25*Math.sqrt(Math.pow(Z,5)/(2*pi*Math.pow(a0,5)))*r*Math.cos(theta)*Math.exp(-Z*r/(2*a0)); } function L_ao_2pz(r,theta,phi,Z) { //Lagrangian of the 2pz atomic orbital return Math.sqrt(Math.pow(Z,7)/(32*pi*Math.pow(a0,7)))*(Z*r/(4*a0)-2)*Math.cos(theta)*Math.exp(-Z*r/(2*a0)); } function ao_3s(r,theta,phi,Z) { //the 3s atomic orbital return 1/81 * Math.sqrt(Math.pow(Z/a0,3) / (3*pi)) * (27 - 18*Z/a0*r + 2*Math.pow(Z/a0*r, 2)) * Math.exp(-Z/a0 * r/3); } function L_ao_3s(r,theta,phi,Z) { //Lagrangian of the 3s atomic orbital return 1/81 * Math.sqrt(Math.pow(Z/a0,3) / (3*pi)) * (-54*Z/a0/r + 39*Math.pow(Z/a0,2) - 6*Math.pow(Z/a0,3)*r + 2/9*Math.pow(Z/a0,4)*r*r) * Math.exp(-Z/a0 * r/3); } function ao_3px(r,theta,phi,Z) { //the 3px atomic orbital return 1/81 * Math.sqrt(2 * Math.pow(Z/a0,3) / pi) * (6 - Z/a0*r) * Z/a0*r * Math.exp(-Z/a0 * r/3) * Math.sin(theta) * Math.cos(phi); } function L_ao_3px(r,theta,phi,Z) { //Lagrangian of the 3px atomic orbital return 1/81 * Math.sqrt(2 * Math.pow(Z/a0,3) / pi) * (-12*Z/a0/r + 8/3*Math.pow(Z/a0,2) - Math.pow(Z/a0,3)*r/9) * Z*r/a0 * Math.exp(-Z/a0 * r/3) * Math.sin(theta) * Math.cos(phi); } function ao_3py(r,theta,phi,Z) { //the 3py atomic orbital return 1/81 * Math.sqrt(2 * Math.pow(Z/a0,3) / pi) * (6-Z/a0*r) * Z/a0*r * Math.exp(-Z/a0 * r/3) * Math.sin(theta) * Math.sin(phi); } function L_ao_3py(r,theta,phi,Z) { //Lagrangian of the 3py atomic orbital return 1/81 * Math.sqrt(2*Math.pow(Z/a0,3) / pi) * (-12*Z/a0/r + 8/3*Math.pow(Z/a0,2) - Math.pow(Z/a0,3)*r/9) * Z*r/a0 * Math.exp(-Z/a0 * r/3) * Math.sin(theta) * Math.sin(phi); } function ao_3pz(r,theta,phi,Z) { //the 3pz atomic orbital return 1/81 * Math.sqrt(2 * Math.pow(Z/a0,3) / pi) * (6-Z/a0*r) * Z/a0*r * Math.exp(-Z/a0 * r/3) * Math.cos(theta); } function L_ao_3pz(r,theta,phi,Z) { //Lagrangian of the 3pz atomic orbital return 1/81 * Math.sqrt(2 * Math.pow(Z/a0,3) / pi) * (-12*Z/a0/r + 8/3*Math.pow(Z/a0,2) - Math.pow(Z/a0,3)*r/9) * Z*r/a0 * Math.exp(-Z/a0 * r/3) * Math.cos(theta); } Z = 1; for (i=0; i<10; i++) { //The energy is calculated for 10 random positions r = Math.random()*a0; theta = Math.random()*pi; phi = Math.random()*2*pi; E = -k1*L_ao_2s(r,theta,phi,Z)/ao_2s(r,theta,phi,Z)-k2*Z/r; document.write(E/ec+" eV<br>"); } </script> (* Content-type: application/vnd.wolfram.mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 10.4' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 158, 7] NotebookDataLength[ 21647, 589] NotebookOptionsPosition[ 20281, 537] NotebookOutlinePosition[ 20639, 553] CellTagsIndexPosition[ 20596, 550] WindowFrame->Normal*) (* Beginning of Notebook Content *) Notebook[{ Cell[CellGroupData[{ Cell[TextData[{ "Atomic orbitals in spherical coordinates\n", StyleBox["Returns atomic orbitals and laplacian of atomic orbitals for \ specified quantum numbers", "Subtitle"] }], "Title", CellChangeTimes->{{3.680501937289125*^9, 3.6805019737784357`*^9}}], Cell[CellGroupData[{ Cell["Evaluate following cell first (shift+enter)", "Chapter", CellChangeTimes->{{3.68053797424739*^9, 3.6805379875702744`*^9}}], Cell[CellGroupData[{ Cell[BoxData[{ RowBox[{ RowBox[{"n", "=", "2"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"l", "=", "1"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"m", "=", "0"}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{"Clear", "[", RowBox[{"a0", ",", "Z"}], "]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n_", ",", " ", "l_", ",", " ", "m_"}], "}"}], ",", " ", "r_", ",", " ", "theta_", ",", "phi_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"2", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "3"}], " ", "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"n", "-", "l", "-", "1"}], "]"}], "/", RowBox[{"(", RowBox[{"2", "*", "n", "*", RowBox[{"Factorial", "[", RowBox[{"n", "+", "l"}], "]"}]}], ")"}]}]}], "]"}], "*", RowBox[{"Exp", "[", RowBox[{ RowBox[{"-", "r"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], "]"}], " ", RowBox[{ RowBox[{"(", RowBox[{"2", "r", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "l"}], " ", "*", RowBox[{"LaguerreL", "[", RowBox[{ RowBox[{"n", " ", "-", " ", "l", " ", "-", " ", "1"}], ",", " ", RowBox[{ RowBox[{"2", " ", "l"}], " ", "+", " ", "1"}], ",", " ", RowBox[{"2", " ", "r", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}]}], "]"}], " ", "*", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"2", "l"}], "+", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"2", "Pi"}], ")"}]}], "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"l", "-", RowBox[{"Abs", "[", "m", "]"}]}], "]"}], "/", RowBox[{"Factorial", "[", RowBox[{"l", "+", RowBox[{"Abs", "[", "m", "]"}]}], "]"}]}]}], "]"}], "*", RowBox[{"LegendreP", "[", RowBox[{"l", ",", " ", RowBox[{"Abs", "[", "m", "]"}], ",", " ", RowBox[{"Cos", "[", "theta", "]"}]}], "]"}], "*", RowBox[{"Sin", "[", RowBox[{ RowBox[{"Abs", "[", "m", "]"}], "*", "phi"}], "]"}]}], "/;", RowBox[{"m", "<", "0"}]}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n_", ",", " ", "l_", ",", " ", "m_"}], "}"}], ",", " ", "r_", ",", " ", "theta_", ",", "phi_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"2", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "3"}], " ", "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"n", "-", "l", "-", "1"}], "]"}], "/", RowBox[{"(", RowBox[{"2", "*", "n", "*", RowBox[{"Factorial", "[", RowBox[{"n", "+", "l"}], "]"}]}], ")"}]}]}], "]"}], "*", RowBox[{"Exp", "[", RowBox[{ RowBox[{"-", "r"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], "]"}], " ", RowBox[{ RowBox[{"(", RowBox[{"2", "r", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "l"}], " ", "*", RowBox[{"LaguerreL", "[", RowBox[{ RowBox[{"n", " ", "-", " ", "l", " ", "-", " ", "1"}], ",", " ", RowBox[{ RowBox[{"2", " ", "l"}], " ", "+", " ", "1"}], ",", " ", RowBox[{"2", " ", "r", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}]}], "]"}], " ", "*", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"2", "l"}], "+", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"4", "Pi"}], ")"}]}], "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"l", "-", "m"}], "]"}], "/", RowBox[{"Factorial", "[", RowBox[{"l", "+", "m"}], "]"}]}]}], "]"}], "*", RowBox[{"LegendreP", "[", RowBox[{"l", ",", " ", "m", ",", " ", RowBox[{"Cos", "[", "theta", "]"}]}], "]"}]}], "/;", RowBox[{"m", "==", "0"}]}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n_", ",", " ", "l_", ",", " ", "m_"}], "}"}], ",", " ", "r_", ",", " ", "theta_", ",", "phi_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"2", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "3"}], " ", "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"n", "-", "l", "-", "1"}], "]"}], "/", RowBox[{"(", RowBox[{"2", "*", "n", "*", RowBox[{"Factorial", "[", RowBox[{"n", "+", "l"}], "]"}]}], ")"}]}]}], "]"}], "*", RowBox[{"Exp", "[", RowBox[{ RowBox[{"-", "r"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], "]"}], " ", RowBox[{ RowBox[{"(", RowBox[{"2", "r", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "l"}], " ", "*", RowBox[{"LaguerreL", "[", RowBox[{ RowBox[{"n", " ", "-", " ", "l", " ", "-", " ", "1"}], ",", " ", RowBox[{ RowBox[{"2", " ", "l"}], " ", "+", " ", "1"}], ",", " ", RowBox[{"2", " ", "r", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}]}], "]"}], " ", "*", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"2", "l"}], "+", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"2", "Pi"}], ")"}]}], "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"l", "-", "m"}], "]"}], "/", RowBox[{"Factorial", "[", RowBox[{"l", "+", "m"}], "]"}]}]}], "]"}], "*", RowBox[{"LegendreP", "[", RowBox[{"l", ",", " ", "m", ",", " ", RowBox[{"Cos", "[", "theta", "]"}]}], "]"}], "*", RowBox[{"Cos", "[", RowBox[{"m", "*", "phi"}], "]"}]}], "/;", RowBox[{"m", ">", "0"}]}]}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"AtOrbit", " ", "=", " ", RowBox[{"FullSimplify", "[", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n", ",", "l", ",", "m"}], "}"}], ",", "r", ",", "\[Theta]", ",", "\[Phi]"}], "]"}], ",", RowBox[{ RowBox[{"Element", "[", RowBox[{ RowBox[{"r", "|", "\[Theta]", "|", "\[Phi]"}], ",", "Reals"}], "]"}], "&&", RowBox[{"r", ">", "0"}], "&&", RowBox[{"\[Theta]", ">", "0"}], "&&", RowBox[{"\[Theta]", "<", "Pi"}], "&&", RowBox[{"\[Phi]", ">", "0"}], "&&", RowBox[{"\[Phi]", "<", RowBox[{"2", "*", "Pi"}]}]}]}], "]"}]}], ";", RowBox[{"LaplaceAtOrbit", " ", "=", " ", RowBox[{"FullSimplify", "[", RowBox[{ RowBox[{"Laplacian", "[", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n", ",", "l", ",", "m"}], "}"}], ",", "r", ",", "\[Theta]", ",", "\[Phi]"}], "]"}], ",", RowBox[{"{", RowBox[{"r", ",", "\[Theta]", ",", "\[Phi]"}], "}"}], ",", "\"\<Spherical\>\""}], "]"}], ",", RowBox[{ RowBox[{"Element", "[", RowBox[{ RowBox[{"r", "|", "\[Theta]", "|", "\[Phi]"}], ",", "Reals"}], "]"}], "&&", RowBox[{"r", ">", "0"}], "&&", RowBox[{"\[Theta]", ">", "0"}], "&&", RowBox[{"\[Theta]", "<", "Pi"}], "&&", RowBox[{"\[Phi]", ">", "0"}], "&&", RowBox[{"\[Phi]", "<", RowBox[{"2", "*", "Pi"}]}]}]}], "]"}]}], ";"}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", "\"\<Atomic orbital\>\"", "\ \[IndentingNewLine]", RowBox[{"AtOrbit", "\[IndentingNewLine]"}], "\[IndentingNewLine]", "\"\<Laplace operator \ applied to atomic orbital displayed above\>\"", "\[IndentingNewLine]", \ "LaplaceAtOrbit"}], "Input", CellChangeTimes->{ 3.35696210375764*^9, 3.3634725675*^9, {3.680267648969605*^9, 3.6802677777893677`*^9}, {3.680267917025701*^9, 3.680267950118266*^9}, { 3.68026803440273*^9, 3.680268078650187*^9}, {3.680268161660758*^9, 3.680268252117462*^9}, {3.680268283719*^9, 3.680268315847728*^9}, { 3.6802683679885*^9, 3.680268401044072*^9}, {3.6802684347322483`*^9, 3.68026846941199*^9}, 3.68026851589177*^9, {3.680268550847946*^9, 3.680268602485694*^9}, {3.680277200598281*^9, 3.680277200967642*^9}, { 3.680277234870984*^9, 3.680277279999036*^9}, {3.680277368558592*^9, 3.680277379167626*^9}, {3.6802775055031853`*^9, 3.6802775433742323`*^9}, { 3.680278118244079*^9, 3.680278127959669*^9}, {3.6802784585585737`*^9, 3.6802785869469976`*^9}, {3.680278865740156*^9, 3.680278911671702*^9}, { 3.6802791705405188`*^9, 3.680279353827753*^9}, {3.6802794897733517`*^9, 3.6802795064406757`*^9}, {3.680279642796474*^9, 3.6802796473713818`*^9}, { 3.68027969209943*^9, 3.680279693816111*^9}, {3.680279868010215*^9, 3.680279944611951*^9}, {3.680280074808546*^9, 3.680280080471993*^9}, { 3.680280136706254*^9, 3.680280136852335*^9}, {3.680280188871863*^9, 3.680280190484304*^9}, {3.680280326138034*^9, 3.680280332936491*^9}, { 3.68028040865331*^9, 3.680280409902307*^9}, {3.680422742828141*^9, 3.680422743874589*^9}, {3.680454127131035*^9, 3.680454139436511*^9}, { 3.6804542844683437`*^9, 3.680454507163357*^9}, {3.68050192461831*^9, 3.680501932191008*^9}, {3.6805019763874397`*^9, 3.680502015775432*^9}, { 3.680502058003353*^9, 3.680502058257691*^9}, {3.6805021194859867`*^9, 3.6805022002999897`*^9}, {3.680538003120837*^9, 3.680538062737301*^9}, { 3.680538113123727*^9, 3.6805381137624083`*^9}}], Cell[BoxData["\<\"Atomic orbital\"\>"], "Output", CellChangeTimes->{{3.680268297943145*^9, 3.680268321703639*^9}, 3.680268371459167*^9, 3.680268402575823*^9, {3.680268571654232*^9, 3.6802686061750917`*^9}, 3.68027709639447*^9, 3.680277202101458*^9, { 3.680277245472929*^9, 3.6802772807869463`*^9}, 3.6802773811394253`*^9, { 3.6802775102931757`*^9, 3.6802775442793913`*^9}, {3.6802781208454733`*^9, 3.680278129004676*^9}, 3.680278590833338*^9, {3.6802788908739767`*^9, 3.680278912749791*^9}, 3.680279171817807*^9, 3.680279355597896*^9, 3.6802795088854847`*^9, 3.6802796481851683`*^9, {3.680279929671982*^9, 3.680279946226309*^9}, 3.680280081912612*^9, 3.680280138253318*^9, 3.680280192226308*^9, {3.680280328045562*^9, 3.680280334796946*^9}, 3.680280411324853*^9, 3.680451203992958*^9, 3.680454484022686*^9, 3.6804545183103657`*^9, {3.680502034951683*^9, 3.6805020601965647`*^9}, 3.680502128398032*^9, {3.680502165483612*^9, 3.680502201131427*^9}, 3.680538087166862*^9, 3.680538117655716*^9}], Cell[BoxData[ FractionBox[ RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"-", FractionBox[ RowBox[{"r", " ", "Z"}], RowBox[{"2", " ", "a0"}]]}]], " ", "r", " ", "Z", " ", SqrtBox[ FractionBox[ SuperscriptBox["Z", "3"], SuperscriptBox["a0", "3"]]], " ", RowBox[{"Cos", "[", "\[Theta]", "]"}]}], RowBox[{"4", " ", "a0", " ", SqrtBox[ RowBox[{"2", " ", "\[Pi]"}]]}]]], "Output", CellChangeTimes->{{3.680268297943145*^9, 3.680268321703639*^9}, 3.680268371459167*^9, 3.680268402575823*^9, {3.680268571654232*^9, 3.6802686061750917`*^9}, 3.68027709639447*^9, 3.680277202101458*^9, { 3.680277245472929*^9, 3.6802772807869463`*^9}, 3.6802773811394253`*^9, { 3.6802775102931757`*^9, 3.6802775442793913`*^9}, {3.6802781208454733`*^9, 3.680278129004676*^9}, 3.680278590833338*^9, {3.6802788908739767`*^9, 3.680278912749791*^9}, 3.680279171817807*^9, 3.680279355597896*^9, 3.6802795088854847`*^9, 3.6802796481851683`*^9, {3.680279929671982*^9, 3.680279946226309*^9}, 3.680280081912612*^9, 3.680280138253318*^9, 3.680280192226308*^9, {3.680280328045562*^9, 3.680280334796946*^9}, 3.680280411324853*^9, 3.680451203992958*^9, 3.680454484022686*^9, 3.6804545183103657`*^9, {3.680502034951683*^9, 3.6805020601965647`*^9}, 3.680502128398032*^9, {3.680502165483612*^9, 3.680502201131427*^9}, 3.680538087166862*^9, 3.680538117656896*^9}], Cell[BoxData["\<\"Laplace operator applied to atomic orbital displayed \ above\"\>"], "Output", CellChangeTimes->{{3.680268297943145*^9, 3.680268321703639*^9}, 3.680268371459167*^9, 3.680268402575823*^9, {3.680268571654232*^9, 3.6802686061750917`*^9}, 3.68027709639447*^9, 3.680277202101458*^9, { 3.680277245472929*^9, 3.6802772807869463`*^9}, 3.6802773811394253`*^9, { 3.6802775102931757`*^9, 3.6802775442793913`*^9}, {3.6802781208454733`*^9, 3.680278129004676*^9}, 3.680278590833338*^9, {3.6802788908739767`*^9, 3.680278912749791*^9}, 3.680279171817807*^9, 3.680279355597896*^9, 3.6802795088854847`*^9, 3.6802796481851683`*^9, {3.680279929671982*^9, 3.680279946226309*^9}, 3.680280081912612*^9, 3.680280138253318*^9, 3.680280192226308*^9, {3.680280328045562*^9, 3.680280334796946*^9}, 3.680280411324853*^9, 3.680451203992958*^9, 3.680454484022686*^9, 3.6804545183103657`*^9, {3.680502034951683*^9, 3.6805020601965647`*^9}, 3.680502128398032*^9, {3.680502165483612*^9, 3.680502201131427*^9}, 3.680538087166862*^9, 3.6805381176602077`*^9}], Cell[BoxData[ FractionBox[ RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"-", FractionBox[ RowBox[{"r", " ", "Z"}], RowBox[{"2", " ", "a0"}]]}]], " ", SuperscriptBox[ RowBox[{"(", FractionBox[ SuperscriptBox["Z", "3"], SuperscriptBox["a0", "3"]], ")"}], RowBox[{"3", "/", "2"}]], " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"-", "8"}], " ", "a0"}], "+", RowBox[{"r", " ", "Z"}]}], ")"}], " ", RowBox[{"Cos", "[", "\[Theta]", "]"}]}], RowBox[{"16", " ", SqrtBox[ RowBox[{"2", " ", "\[Pi]"}]], " ", "Z"}]]], "Output", CellChangeTimes->{{3.680268297943145*^9, 3.680268321703639*^9}, 3.680268371459167*^9, 3.680268402575823*^9, {3.680268571654232*^9, 3.6802686061750917`*^9}, 3.68027709639447*^9, 3.680277202101458*^9, { 3.680277245472929*^9, 3.6802772807869463`*^9}, 3.6802773811394253`*^9, { 3.6802775102931757`*^9, 3.6802775442793913`*^9}, {3.6802781208454733`*^9, 3.680278129004676*^9}, 3.680278590833338*^9, {3.6802788908739767`*^9, 3.680278912749791*^9}, 3.680279171817807*^9, 3.680279355597896*^9, 3.6802795088854847`*^9, 3.6802796481851683`*^9, {3.680279929671982*^9, 3.680279946226309*^9}, 3.680280081912612*^9, 3.680280138253318*^9, 3.680280192226308*^9, {3.680280328045562*^9, 3.680280334796946*^9}, 3.680280411324853*^9, 3.680451203992958*^9, 3.680454484022686*^9, 3.6804545183103657`*^9, {3.680502034951683*^9, 3.6805020601965647`*^9}, 3.680502128398032*^9, {3.680502165483612*^9, 3.680502201131427*^9}, 3.680538087166862*^9, 3.6805381176615143`*^9}] }, Open ]], Cell[CellGroupData[{ Cell[TextData[{ StyleBox["Evaluate following cells to get ", "Chapter"], Cell[BoxData[ FormBox[ StyleBox[ RowBox[{"L", StyleBox[ AdjustmentBox["A", BoxBaselineShift->-0.4, BoxMargins->{{-0.5, -0.3}, {0, 0}}], FontSize->Smaller], "T", AdjustmentBox["E", BoxBaselineShift->0.5, BoxMargins->{{-0.3, 0}, {0, 0}}], "X"}], SingleLetterItalics->False], TraditionalForm]]], StyleBox[" syntax", "Chapter"], " " }], "Section", CellChangeTimes->{{3.680537047447669*^9, 3.6805370920176907`*^9}, { 3.68053741912608*^9, 3.680537449967037*^9}, {3.680537552860992*^9, 3.680537560032627*^9}}], Cell[CellGroupData[{ Cell[BoxData[{"\"\<Atomic Orbital\>\"", "\[IndentingNewLine]", RowBox[{"TeXForm", "[", "AtOrbit", "]"}]}], "Input", CellChangeTimes->{{3.680537452108548*^9, 3.680537509094211*^9}, { 3.680537548839533*^9, 3.6805375672542763`*^9}}], Cell[BoxData["\<\"Atomic Orbital\"\>"], "Output", CellChangeTimes->{ 3.6805375694238787`*^9, {3.6805380947972107`*^9, 3.680538121837367*^9}}], Cell["\<\ \\frac{r Z \\sqrt{\\frac{Z^3}{\\text{a0}^3}} \\cos (\\theta ) e^{-\\frac{r Z}{2 \\text{a0}}}}{4 \\sqrt{2 \\pi } \\text{a0}}\ \>", "Output", CellChangeTimes->{ 3.6805375694238787`*^9, {3.6805380947972107`*^9, 3.6805381218384666`*^9}}] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{"\[IndentingNewLine]", RowBox[{ "\"\<Laplace operator applied to atomic orbital displayed above\>\"", "\[IndentingNewLine]", RowBox[{"TeXForm", "[", "LaplaceAtOrbit", "]"}]}]}]], "Input", CellChangeTimes->{{3.680537452108548*^9, 3.680537509094211*^9}, { 3.680537548839533*^9, 3.6805375672542763`*^9}}], Cell[BoxData["\<\"Laplace operator applied to atomic orbital displayed \ above\"\>"], "Output", CellChangeTimes->{3.680537511655724*^9, 3.6805375723269787`*^9, 3.68053812449619*^9}], Cell["\<\ \\frac{\\left(\\frac{Z^3}{\\text{a0}^3}\\right)^{3/2} \\cos (\\theta ) e^{-\\frac{r Z}{2 \\text{a0}}} (r Z-8 \\text{a0})}{16 \\sqrt{2 \\pi } Z}\ \>", "Output", CellChangeTimes->{3.680537511655724*^9, 3.6805375723269787`*^9, 3.680538124497911*^9}] }, Open ]] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["\<\ Evaluate following cells to get orbitals in mathematica syntax\ \>", "Chapter", CellChangeTimes->{{3.6805375953984947`*^9, 3.680537709515551*^9}}], Cell[CellGroupData[{ Cell[BoxData[{"\"\<Atomic Orbital\>\"", "\[IndentingNewLine]", RowBox[{"InputForm", "[", "AtOrbit", "]"}]}], "Input", CellChangeTimes->{{3.6805375953984947`*^9, 3.6805376593171663`*^9}, 3.6805377145532017`*^9}], Cell[BoxData["\<\"Atomic Orbital\"\>"], "Output", CellChangeTimes->{3.6805377176185207`*^9, 3.680538128774366*^9}], Cell["\<\ (r*Z*Sqrt[Z^3/a0^3]*Cos[\[Theta]])/(4*a0*E^((r*Z)/(2*a0))* Sqrt[2*Pi])\ \>", "Output", CellChangeTimes->{3.6805377176185207`*^9, 3.680538128776547*^9}] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[{"\"\<Laplace operator applied to atomic orbital displayed above\ \>\"", "\[IndentingNewLine]", RowBox[{"InputForm", "[", "LaplaceAtOrbit", "]"}]}], "Input", CellChangeTimes->{{3.6805375953984947`*^9, 3.6805376593171663`*^9}, 3.6805377145532017`*^9}], Cell[BoxData["\<\"Laplace operator applied to atomic orbital displayed \ above\"\>"], "Output", CellChangeTimes->{3.680537719968749*^9, 3.680538131919155*^9}], Cell["\<\ ((Z^3/a0^3)^(3/2)*(-8*a0 + r*Z)*Cos[\[Theta]])/ (16*E^((r*Z)/(2*a0))*Sqrt[2*Pi]*Z)\ \>", "Output", CellChangeTimes->{3.680537719968749*^9, 3.680538131921568*^9}] }, Open ]] }, Open ]] }, Open ]] }, WindowSize->{1280, 1376}, WindowMargins->{{Automatic, -10}, {Automatic, -8}}, Magnification->1.5, FrontEndVersion->"10.4 for Linux x86 (64-bit) (April 11, 2016)", StyleDefinitions->"Default.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[580, 22, 257, 5, 239, "Title"], Cell[CellGroupData[{ Cell[862, 31, 129, 1, 98, "Chapter"], Cell[CellGroupData[{ Cell[1016, 36, 10456, 280, 1039, "Input"], Cell[11475, 318, 1042, 14, 45, "Output"], Cell[12520, 334, 1437, 29, 131, "Output"], Cell[13960, 365, 1090, 15, 45, "Output"], Cell[15053, 382, 1618, 36, 106, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[16708, 423, 648, 20, 109, "Section"], Cell[CellGroupData[{ Cell[17381, 447, 234, 3, 75, "Input"], Cell[17618, 452, 144, 2, 45, "Output"], Cell[17765, 456, 252, 6, 111, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[18054, 467, 340, 7, 104, "Input"], Cell[18397, 476, 185, 3, 45, "Output"], Cell[18585, 481, 266, 6, 111, "Output"] }, Open ]] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[18912, 494, 157, 3, 165, "Chapter"], Cell[CellGroupData[{ Cell[19094, 501, 217, 3, 75, "Input"], Cell[19314, 506, 115, 1, 78, "Output"], Cell[19432, 509, 164, 4, 87, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[19633, 518, 270, 4, 75, "Input"], Cell[19906, 524, 159, 2, 45, "Output"], Cell[20068, 528, 173, 4, 87, "Output"] }, Open ]] }, Open ]] }, Open ]] } ] *) #include <iostream> #include <cmath> #include <time.h> using namespace std; double phi1s(double r,double t,double p); double Lphi1s(double r,double t,double p); double phi2s(double r,double t,double p); double Lphi2s(double r,double t,double p); double phi2px(double r,double t,double p); double Lphi2px(double r,double t,double p); double phi2py(double r,double t,double p); double Lphi2py(double r,double t,double p); double phi2pz(double r,double t,double p); double Lphi2pz(double r,double t,double p); double phi3s(double r,double t,double p); double Lphi3s(double r,double t,double p); double phi3px(double r,double t,double p); double Lphi3px(double r,double t,double p); double phi3py(double r,double t,double p); double Lphi3py(double r,double t,double p); double phi3pz(double r,double t,double p); double Lphi3pz(double r,double t,double p); double Z=1; double pi=3.141592; int main() { double r,t,p; double x1=0,y1=10,x2=0,y2=pi,x3=0,y3=2*pi; srand (time(NULL)); cout << "This program tests if the hydrogen-wavefunctions evaluated on random positions (r,t,p) in spherical coordinates lead to the correct energies E=H|phi_n(r,t,p)/phi_n(r,t,p)=-Z^2/n^2" << endl << endl; cout << "The proram uses atomic units where the energy is given in 1 Hartree = 27.21 eV, = 2 Rydberg " << endl << endl; cout << "E / Hartree : " << endl << endl; cout << "phi1s(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi1s(r,t,p)-phi1s(r,t,p)*Z/r)/phi1s(r,t,p) << endl; } cout << endl; cout << "phi2s(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi2s(r,t,p)-phi2s(r,t,p)*Z/r)/phi2s(r,t,p) << endl; } cout << endl; cout << "phi2px(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi2px(r,t,p)-phi2px(r,t,p)*Z/r)/phi2px(r,t,p) << endl; } cout << endl; cout << "phi2py(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi2py(r,t,p)-phi2py(r,t,p)*Z/r)/phi2py(r,t,p) << endl; } cout << endl; cout << "phi2pz(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi2pz(r,t,p)-phi2pz(r,t,p)*Z/r)/phi2pz(r,t,p) << endl; } cout << endl; cout << "phi3s(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi3s(r,t,p)-phi3s(r,t,p)*Z/r)/phi3s(r,t,p) << endl; } cout << endl; cout << "phi3px(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi3px(r,t,p)-phi3px(r,t,p)*Z/r)/phi3px(r,t,p) << endl; } cout << endl; cout << "phi3py(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi3py(r,t,p)-phi3py(r,t,p)*Z/r)/phi3py(r,t,p) << endl; } cout << endl; cout << "phi3pz(r,t,p): " << endl; for(int i=0; i<10; i++) { r=((rand()%100)/100.0)*(y1-x1)+1.0E-10; t=((rand()%100)/100.0)*(y2-x2)+1.0E-10; p=((rand()%100)/100.0)*(y3-x3)+1.0E-10; cout << "r=" << r << ", " << "t=" << t << ", " << "p=" << p << " : "; cout << (-0.5*Lphi3pz(r,t,p)-phi3pz(r,t,p)*Z/r)/phi3pz(r,t,p) << endl; } cout << endl; int wait; cin >> wait; } double phi1s(double r,double t,double p) { return sqrt(pow(Z,3)/pi)*exp(-Z*r); } double Lphi1s(double r,double t,double p) { return sqrt(pow(Z,5)/pi)*(Z-2.0/r)*exp(-Z*r); } double phi2s(double r,double t,double p) { return (2.0-r)*exp(-r/2.0)/(4.0*sqrt(2.0*pi)); } double Lphi2s(double r,double t,double p) { return (10.0*r-16.0-r*r)*exp(-r/2.0)/(16.0*r*sqrt(2.0*pi)); } double phi2px(double r,double t,double p) { return -sqrt(pow(Z,3)/pi)*Z*r*exp(-Z*r/2.0)*sin(t)*cos(p)/8.0; } double Lphi2px(double r,double t,double p) { return sqrt(pow(Z,7)/pi)*(8.0-Z*r)*exp(-Z*r/2.0)*cos(p)*sin(t)/32.0; } double phi2py(double r,double t,double p) { return sqrt(pow(Z,3)/pi)*Z*r*exp(-Z*r/2.0)*sin(t)*sin(p)/8.0; } double Lphi2py(double r,double t,double p) { return sqrt(pow(Z,7)/pi)*(Z*r-8.0)*exp(-Z*r/2.0)*sin(p)*sin(t)/32.0; } double phi2pz(double r,double t,double p) { return sqrt(pow(Z,3)/(2.0*pi))*Z*r*exp(-Z*r/2.0)*cos(t)/4.0; } double Lphi2pz(double r,double t,double p) { return sqrt(pow(Z,7)/(2.0*pi))*(Z*r-8)*exp(-Z*r/2.0)*cos(t)/16.0; } double phi3s(double r,double t,double p) { return sqrt(pow(Z,3)/(3.0*pi))*(27.0-18.0*Z*r+2.0*Z*Z*r*r)*exp(-Z*r/3.0)/81.0; } double Lphi3s(double r,double t,double p) { return sqrt(pow(Z,5)/(3.0*pi))*(2.0*r*r*r*Z*Z*Z-54.0*r*r*Z*Z+351.0*r*Z-486.0)*exp(-Z*r/3.0)/(729.0*r); } double phi3px(double r,double t,double p) { return sqrt(2.0*pow(Z,3)/pi)*Z*r*(6.0-Z*r)*exp(-Z*r/3.0)*sin(t)*cos(p)/81.0; } double Lphi3px(double r,double t,double p) { return sqrt(2.0*pow(Z,7)/pi)*(6.0-Z*r)*(Z*r-18.0)*exp(-Z*r/3.0)*cos(p)*sin(t)/729.0; } double phi3py(double r,double t,double p) { return sqrt(2.0*pow(Z,3)/pi)*Z*r*(6.0-Z*r)*exp(-Z*r/3.0)*sin(t)*sin(p)/81.0; } double Lphi3py(double r,double t,double p) { return -sqrt(2.0*pow(Z,7)/pi)*(108.0-24.0*Z*r+Z*Z*r*r)*exp(-Z*r/3.0)*sin(t)*sin(p)/729.0; } double phi3pz(double r,double t,double p) { return sqrt(2.0*pow(Z,3)/pi)*Z*r*(6.0-Z*r)*exp(-Z*r/3.0)*cos(t)/81.0; } double Lphi3pz(double r,double t,double p) { return -sqrt(2.0*pow(Z,7)/pi)*(108.0-24.0*Z*r+Z*Z*r*r)*exp(-Z*r/3.0)*cos(t)/729.0; }

If an atom is located at position $(x_i,y_i,z_i)$, sometimes it is easier to work in Cartesan coordinates where $r=\sqrt{(x-x_i)^2+(y-y_i)^2+(z-z_i)^2}$.

\[\begin{equation} \begin{array}{c,c} \phi_{1s}=\sqrt{\frac{Z^3}{\pi a^3_0}}\exp\left(-\frac{Zr}{a_0}\right)& \nabla^2\phi_{1s}=\sqrt{\frac{Z^5}{\pi a^5_0}} \left(\frac{Z}{a_0}-\frac{2}{r}\right)\exp\left(-\frac{Zr}{a_0}\right)\\ \phi_{2s}=\frac{1}{4} \sqrt{\frac{Z^3}{2\pi a_0^3}}\left(2-\frac{Zr}{a_0}\right)\exp\left(-\frac{Zr}{2a_0}\right)& \nabla^2\phi_{2s}=\frac{1}{4}\sqrt{\frac{Z^5}{2\pi a^5_0}} \left(-\frac{Z^2r}{4a_0^2}+\frac{5Z}{2a_0}-\frac{4}{r}\right)\exp\left(-\frac{Zr}{2a_0}\right)\\ \phi_{2p_x}=\frac{1}{4} \sqrt{\frac{Z^5}{2\pi a_0^5}}(x-x_i) \exp\left(-\frac{Zr}{2a_0}\right) & \nabla^2\phi_{2p_x}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\frac{x-x_i}{r} \exp\left(-\frac{Zr}{2a_0}\right),\\ \phi_{2p_y}=\frac{1}{4} \sqrt{\frac{Z^5}{2\pi a_0^5}}(y-y_i) \exp\left(-\frac{Zr}{2a_0}\right) & \nabla^2\phi_{2p_y}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\frac{y-y_i}{r} \exp\left(-\frac{Zr}{2a_0}\right),\\ \phi_{2p_z}=\frac{1}{4} \sqrt{\frac{Z^5}{2\pi a_0^5}}(z-z_i) \exp\left(-\frac{Zr}{2a_0}\right) & \nabla^2\phi_{2p_z}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\frac{z-z_i}{r} \exp\left(-\frac{Zr}{2a_0}\right),\\ \phi_{3s} = \frac{1}{81}\sqrt{\frac{Z^3}{3\pi a_0^3}} \left(27-\frac{18Zr}{a_0}+\frac{2Z^2r^2}{a_0^2}\right) \exp\left(\frac{-Zr}{3a_0}\right) & \nabla^2\phi_{3s} = \frac{1}{81}\sqrt{\frac{Z^3}{3\pi a_0^3}} \left(-\frac{54 Z}{a_0 r}+\frac{39Z^2}{a_0^2}-\frac{6Z^3r}{a_0^3}+\frac{2Z^4r^2}{9a_0^4}\right) \exp\left(\frac{-Zr}{3a_0}\right) \\ \phi_{3p_x} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(6-\frac{Zr}{a_0}\right)\frac{Z(x-x_i)}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) & \nabla^2\phi_{3p_x} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(-\frac{12Z}{a_0 r}+\frac{8Z^2}{3a_0^2}-\frac{Z^3 r}{9a_0^3}\right) \frac{Z(x-x_i)}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \\ \phi_{3p_y} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(6-\frac{Zr}{a_0}\right)\frac{Z(y-y_i)}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) & \nabla^2\phi_{3p_y} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(-\frac{12Z}{a_0 r}+\frac{8Z^2}{3a_0^2}-\frac{Z^3r}{9a_0^3}\right) \frac{Z(y-y_i)}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \\ \phi_{3p_z} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(6-\frac{Zr}{a_0}\right)\frac{Z(z-z_i)}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) & \nabla^2\phi_{3p_z} = \frac{1}{81}\sqrt{\frac{2Z^3}{\pi a_0^3}} \left(-\frac{12Z}{a_0 r}+\frac{8Z^2}{3a_0^2}-\frac{Z^3 r}{9a_0^3}\right) \frac{Z(z-z_i)}{a_0} \exp\left(\frac{-Zr}{3a_0}\right) \\ \phi_{3d_{z^2}} = \frac{1}{81}\sqrt{\frac{Z^7}{6\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\left(3\frac{(z-z_i)^2}{r^2}-1\right) & \nabla^2\phi_{3d_{z^2}} = \frac{1}{243}\sqrt{\frac{Z^9}{6\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\left(3\frac{(z-z_i)^2}{r^2}-1\right)\\ \phi_{3d_{xz}} = \frac{1}{81}\sqrt{\frac{2Z^7}{\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\frac{(x-x_i)(z-z_i)}{r^2} & \nabla^2\phi_{3d_{xz}} = \frac{1}{243}\sqrt{\frac{2Z^9}{\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\frac{(x-x_i)(z-z_i)}{r^2}\\ \phi_{3d_{yz}} = \frac{1}{81}\sqrt{\frac{2Z^7}{\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\frac{(y-y_i)(z-z_i)}{r^2} & \nabla^2\phi_{3d_{yz}} = \frac{1}{243}\sqrt{\frac{2Z^9}{\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\frac{(y-y_i)(z-z_i)}{r^2}\\ \phi_{3d_{xy}} = \frac{1}{81}\sqrt{\frac{2Z^7}{\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\frac{(x-x_i)(y-y_i)}{r^2} & \nabla^2\phi_{3d_{xy}} = \frac{1}{243}\sqrt{\frac{2Z^9}{\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\frac{(x-x_i)(y-y_i)}{r^2}\\ \phi_{3d_{x^2-y^2}} = \frac{1}{81}\sqrt{\frac{Z^7}{2\pi a_0^7}}~r^2\exp\left(\frac{-Zr}{3a_0}\right)\frac{(x-x_i)^2-(y-y_i)^2}{r^2} & \nabla^2\phi_{3d_{x^2-y^2}} = \frac{1}{243}\sqrt{\frac{Z^9}{2\pi a_0^9}}~r\left(\frac{Z}{3a_0}r-6\right)\exp\left(\frac{-Zr}{3a_0}\right)\frac{(x-x_i)^2-(y-y_i)^2}{r^2}\\ \end{array} \end{equation}\]

The code below defines the atomic orbitals in Cartesian coordinates centered at position $(x_i,y_i,z_i)$.

<script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = 3.1415926535897932; var hbar = 1.05457266E-34; var ec = 1.60217733E-19; var eps0 = 8.854187817E-12; var me = 9.1093897E-31; var k1 = 0.5*hbar*hbar/me; var k2 = ec*ec/(4*pi*eps0); //Coulomb constant function ao_1s(x,y,z,xi,yi,zi,Za) { //the 1s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.sqrt(Math.pow(Za,3)/(pi*Math.pow(a0,3)))*Math.exp(-Za*r/a0); } function L_ao_1s(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 1s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.sqrt(Math.pow(Za,5)/(pi*Math.pow(a0,5)))*(Za/a0 -2/r)*Math.exp(-Za*r/a0); } function ao_2s(x,y,z,xi,yi,zi,Za) { //the 2s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return 0.25*Math.sqrt(Math.pow(Za,3)/(2*pi*Math.pow(a0,3)))*(2-Za*r/a0)*Math.exp(-Za*r/(2*a0)); } function L_ao_2s(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 2s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return 0.25*Math.sqrt(Math.pow(Za,5)/(2*pi*Math.pow(a0,5)))*(-Za*Za*r/(4*a0*a0)+5*Za/(2*a0)-4/r)*Math.exp(-Za*r/(2*a0)); } function ao_2pz(x,y,z,xi,yi,zi,Za) { //the 2pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.sqrt(Math.pow(Za/(2*a0),5)/pi)*(z-zi)*Math.exp(-Za*r/(2*a0)); } function L_ao_2pz(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 2pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.pow(Za/a0,7/2)*Math.pow(1/2,5/2)*Math.sqrt(1/pi)*(Za*r/(4*a0)-2)*((z-zi)/r)*Math.exp(-Za*r/(2*a0)); } function ao_3s(x,y,z,xi,yi,zi,Za) { //the 3s atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(Math.pow(Za/a0,3) / (3*pi)) * (27 - 18*Za/a0*r + 2*Math.pow(Za/a0*r, 2)) * Math.exp(-Za/a0 * r/3); } function L_ao_3s(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3s atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(Math.pow(Za/a0,3) / (3*pi)) * (-54*Za/a0/r + 39*Math.pow(Za/a0,2) - 6*Math.pow(Za/a0,3)*r + 2/9*Math.pow(Za/a0,4)*r*r) * Math.exp(-Za/a0 * r/3); } function ao_3px(x,y,z,xi,yi,zi,Za) { //the 3px atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (6 - Za/a0*r) * Za/a0*(x-xi) * Math.exp(-Za/a0 * r/3); } function L_ao_3px(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3px atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (-12*Za/a0/r + 8/3*Math.pow(Za/a0,2) - Math.pow(Za/a0,3)*r/9) * Za/a0*(x-xi) * Math.exp(-Za/a0 * r/3); } function ao_3py(x,y,z,xi,yi,zi,Za) { //the 3py atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (6 - Za/a0*r) * Za/a0*(y-yi) * Math.exp(-Za/a0 * r/3); } function L_ao_3py(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3py atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (-12*Za/a0/r + 8/3*Math.pow(Za/a0,2) - Math.pow(Za/a0,3)*r/9) * Za/a0*(y-yi) * Math.exp(-Za/a0 * r/3); } function ao_3pz(x,y,z,xi,yi,zi,Za) { //the 3pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (6 - Za/a0*r) * Za/a0*(z-zi) * Math.exp(-Za/a0 * r/3); } function L_ao_3pz(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (-12*Za/a0/r + 8/3*Math.pow(Za/a0,2) - Math.pow(Za/a0,3)*r/9) * Za/a0*(z-zi) * Math.exp(-Za/a0 * r/3); } function ao_3dxz(x,y,z,xi,yi,zi,Za) { //the 3dxz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2*Za**7/(pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((x-xi)*(z-zi))/r**2; } function L_ao_3dxz(x,y,z,xi,yi,zi,Za) { //the 3dxz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(2*Za**9/(pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((x-xi)*(z-zi))/r**2; } function ao_3dxy(x,y,z,xi,yi,zi,Za) { //the 3dxy atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2*Za**7/(pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((x-xi)*(y-yi))/r**2; } function L_ao_3dxy(x,y,z,xi,yi,zi,Za) { //the 3dxy atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(2*Za**9/(pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((x-xi)*(y-yi))/r**2; } function ao_3dyz(x,y,z,xi,yi,zi,Za) { //the 3dyz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2*Za**7/(pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((y-yi)*(z-zi))/r**2; } function L_ao_3dyz(x,y,z,xi,yi,zi,Za) { //the 3dyz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(2*Za**9/(pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((y-yi)*(z-zi))/r**2; } function ao_3dzz(x,y,z,xi,yi,zi,Za) { //the 3dz^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81*Math.sqrt(Za**7/(6*pi*a0**7))*Math.exp(-Za*r/(3*a0))*r**2*(3*(z-zi)**2/r**2-1); } function L_ao_3dzz(x,y,z,xi,yi,zi,Za) { //the 3dz^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243*Math.sqrt(Za**9/(6*pi*a0**9))*Math.exp(-Za*r/(3*a0))*r*(Za*r/(3*a0)-6)*(3*(z-zi)**2/r**2-1); } function ao_3dx2my2(x,y,z,xi,yi,zi,Za) { //the 3dx^2-y^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(Za**7/(2*pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((x-xi)*(x-xi)-(y-yi)*(y-yi))/r**2; } function L_ao_3dx2my2(x,y,z,xi,yi,zi,Za) { //the 3dx^2-y^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(Za**9/(2*pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((x-xi)*(x-xi)-(y-yi)*(y-yi))/r**2; } Za = 1; for (i=0; i<10; i++) { //The energy is calculated for 10 random positions x = 2*a0*(0.5-Math.random()); xi = 2*a0*(0.5-Math.random()); y = 2*a0*(0.5-Math.random()); yi = 2*a0*(0.5-Math.random()); z = 2*a0*(0.5-Math.random()); zi = 2*a0*(0.5-Math.random()); E = -k1*L_ao_2s(x,y,z,xi,yi,zi,Za)/ao_2s(x,y,z,xi,yi,zi,Za)-k2*Za/r; document.write(E/ec+" eV<br>"); } </script> (* Content-type: application/vnd.wolfram.mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 10.4' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 158, 7] NotebookDataLength[ 33186, 851] NotebookOptionsPosition[ 31666, 794] NotebookOutlinePosition[ 32023, 810] CellTagsIndexPosition[ 31980, 807] WindowFrame->Normal*) (* Beginning of Notebook Content *) Notebook[{ Cell[BoxData[{ StyleBox[ StyleBox[ RowBox[{ "Atomic", " ", "orbitals", " ", "in", " ", "cartesian", " ", "coordinates"}], "Title"], FontWeight->"Plain"], "\[IndentingNewLine]", StyleBox[ StyleBox[ RowBox[{ "Returns", " ", "atomic", " ", "orbitals", " ", "and", " ", "laplacian", " ", "of", " ", "atomic", " ", "orbitals", " ", "for", " ", "specified", " ", "quantum", " ", "numbers"}], "Subtitle"], FontWeight->"Plain"]}], "Input", CellChangeTimes->{ 3.35696210375764*^9, 3.3634725675*^9, {3.680267648969605*^9, 3.6802677777893677`*^9}, {3.680267917025701*^9, 3.680267950118266*^9}, { 3.68026803440273*^9, 3.680268078650187*^9}, {3.680268161660758*^9, 3.680268252117462*^9}, {3.680268283719*^9, 3.680268315847728*^9}, { 3.6802683679885*^9, 3.680268401044072*^9}, {3.6802684347322483`*^9, 3.68026846941199*^9}, 3.68026851589177*^9, {3.680268550847946*^9, 3.680268602485694*^9}, {3.680277200598281*^9, 3.680277200967642*^9}, { 3.680277234870984*^9, 3.680277279999036*^9}, {3.680277368558592*^9, 3.680277379167626*^9}, {3.6802775055031853`*^9, 3.6802775433742323`*^9}, { 3.680278118244079*^9, 3.680278127959669*^9}, {3.6802784585585737`*^9, 3.6802785869469976`*^9}, {3.680278865740156*^9, 3.680278911671702*^9}, { 3.6802791705405188`*^9, 3.680279353827753*^9}, {3.6802794897733517`*^9, 3.6802795064406757`*^9}, {3.680279642796474*^9, 3.6802796473713818`*^9}, { 3.68027969209943*^9, 3.680279693816111*^9}, {3.680279868010215*^9, 3.680279944611951*^9}, {3.680280074808546*^9, 3.680280080471993*^9}, { 3.680280136706254*^9, 3.680280136852335*^9}, {3.680280188871863*^9, 3.680280190484304*^9}, {3.680280326138034*^9, 3.680280332936491*^9}, { 3.68028040865331*^9, 3.680280409902307*^9}, {3.680422742828141*^9, 3.680422743874589*^9}, {3.680454127131035*^9, 3.680454139436511*^9}, { 3.6804542844683437`*^9, 3.680454507163357*^9}, {3.6804982821582613`*^9, 3.6804983117917747`*^9}, {3.680498545809291*^9, 3.6804986280861063`*^9}, { 3.680498859127798*^9, 3.6804990405705233`*^9}, {3.680499458218275*^9, 3.6804995193350782`*^9}, {3.680499994141029*^9, 3.680500114217901*^9}, { 3.680500145428303*^9, 3.6805001801054068`*^9}, {3.680500356902849*^9, 3.680500381139038*^9}, {3.6805004478266573`*^9, 3.68050070843865*^9}, 3.68050080824401*^9, {3.6805371041846943`*^9, 3.680537146537182*^9}}], Cell[CellGroupData[{ Cell["Evaluate following cell first (shift + enter)", "Chapter", CellChangeTimes->{ 3.35696210375764*^9, 3.3634725675*^9, {3.680267648969605*^9, 3.6802677777893677`*^9}, {3.680267917025701*^9, 3.680267950118266*^9}, { 3.68026803440273*^9, 3.680268078650187*^9}, {3.680268161660758*^9, 3.680268252117462*^9}, {3.680268283719*^9, 3.680268315847728*^9}, { 3.6802683679885*^9, 3.680268401044072*^9}, {3.6802684347322483`*^9, 3.68026846941199*^9}, 3.68026851589177*^9, {3.680268550847946*^9, 3.680268602485694*^9}, {3.680277200598281*^9, 3.680277200967642*^9}, { 3.680277234870984*^9, 3.680277279999036*^9}, {3.680277368558592*^9, 3.680277379167626*^9}, {3.6802775055031853`*^9, 3.6802775433742323`*^9}, { 3.680278118244079*^9, 3.680278127959669*^9}, {3.6802784585585737`*^9, 3.6802785869469976`*^9}, {3.680278865740156*^9, 3.680278911671702*^9}, { 3.6802791705405188`*^9, 3.680279353827753*^9}, {3.6802794897733517`*^9, 3.6802795064406757`*^9}, {3.680279642796474*^9, 3.6802796473713818`*^9}, { 3.68027969209943*^9, 3.680279693816111*^9}, {3.680279868010215*^9, 3.680279944611951*^9}, {3.680280074808546*^9, 3.680280080471993*^9}, { 3.680280136706254*^9, 3.680280136852335*^9}, {3.680280188871863*^9, 3.680280190484304*^9}, {3.680280326138034*^9, 3.680280332936491*^9}, { 3.68028040865331*^9, 3.680280409902307*^9}, {3.680422742828141*^9, 3.680422743874589*^9}, {3.680454127131035*^9, 3.680454139436511*^9}, { 3.6804542844683437`*^9, 3.680454507163357*^9}, {3.6804982821582613`*^9, 3.6804983117917747`*^9}, {3.680498545809291*^9, 3.6804986280861063`*^9}, { 3.680498859127798*^9, 3.6804990405705233`*^9}, {3.680499458218275*^9, 3.6804995193350782`*^9}, {3.680499994141029*^9, 3.680500114217901*^9}, { 3.680500145428303*^9, 3.6805001801054068`*^9}, {3.680500356902849*^9, 3.680500381139038*^9}, {3.6805004478266573`*^9, 3.68050070843865*^9}, 3.68050080824401*^9, {3.6805371041846943`*^9, 3.6805371921655407`*^9}}], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{"\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"n", "=", "2"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"l", "=", "1"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"m", "=", "0"}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"Clear", "[", RowBox[{"a0", ",", "Z", ",", "r"}], "]"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n_", ",", " ", "l_", ",", " ", "m_"}], "}"}], ",", " ", "x_", ",", " ", "y_", ",", "z_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"2", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "3"}], " ", "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"n", "-", "l", "-", "1"}], "]"}], "/", RowBox[{"(", RowBox[{"2", "*", "n", "*", RowBox[{"Factorial", "[", RowBox[{"n", "+", "l"}], "]"}]}], ")"}]}]}], "]"}], "*", RowBox[{"Exp", "[", RowBox[{ RowBox[{"-", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}]}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], "]"}], " ", RowBox[{ RowBox[{"(", RowBox[{"2", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "l"}], " ", "*", RowBox[{"LaguerreL", "[", RowBox[{ RowBox[{"n", " ", "-", " ", "l", " ", "-", " ", "1"}], ",", " ", RowBox[{ RowBox[{"2", " ", "l"}], " ", "+", " ", "1"}], ",", " ", RowBox[{"2", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}]}], "]"}], " ", "*", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"2", "l"}], "+", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"2", "Pi"}], ")"}]}], "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"l", "-", RowBox[{"Abs", "[", "m", "]"}]}], "]"}], "/", RowBox[{"Factorial", "[", RowBox[{"l", "+", RowBox[{"Abs", "[", "m", "]"}]}], "]"}]}]}], "]"}], "*", RowBox[{"LegendreP", "[", RowBox[{"l", ",", " ", RowBox[{"Abs", "[", "m", "]"}], ",", " ", RowBox[{"z", "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}]}]}], "]"}], "*", RowBox[{"Sin", "[", RowBox[{ RowBox[{"Abs", "[", "m", "]"}], "*", RowBox[{"ArcSin", "[", RowBox[{"y", "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}]}], "]"}]}], "]"}]}], "]"}]}], "/;", RowBox[{"m", "<", "0"}]}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n_", ",", " ", "l_", ",", " ", "m_"}], "}"}], ",", " ", "x_", ",", " ", "y_", ",", "z_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"2", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "3"}], " ", "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"n", "-", "l", "-", "1"}], "]"}], "/", RowBox[{"(", RowBox[{"2", "*", "n", "*", RowBox[{"Factorial", "[", RowBox[{"n", "+", "l"}], "]"}]}], ")"}]}]}], "]"}], "*", RowBox[{"Exp", "[", RowBox[{ RowBox[{"-", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}]}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], "]"}], " ", RowBox[{ RowBox[{"(", RowBox[{"2", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "l"}], " ", "*", RowBox[{"LaguerreL", "[", RowBox[{ RowBox[{"n", " ", "-", " ", "l", " ", "-", " ", "1"}], ",", " ", RowBox[{ RowBox[{"2", " ", "l"}], " ", "+", " ", "1"}], ",", " ", RowBox[{"2", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}]}], "]"}], " ", "*", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"2", "l"}], "+", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"4", "Pi"}], ")"}]}], "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"l", "-", "m"}], "]"}], "/", RowBox[{"Factorial", "[", RowBox[{"l", "+", "m"}], "]"}]}]}], "]"}], "*", RowBox[{"LegendreP", "[", RowBox[{"l", ",", " ", RowBox[{"Abs", "[", "m", "]"}], ",", " ", RowBox[{"z", "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}]}]}], "]"}]}], "/;", RowBox[{"m", "==", "0"}]}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n_", ",", " ", "l_", ",", " ", "m_"}], "}"}], ",", " ", "x_", ",", " ", "y_", ",", "z_"}], "]"}], " ", ":=", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{"2", "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "3"}], " ", "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"n", "-", "l", "-", "1"}], "]"}], "/", RowBox[{"(", RowBox[{"2", "*", "n", "*", RowBox[{"Factorial", "[", RowBox[{"n", "+", "l"}], "]"}]}], ")"}]}]}], "]"}], "*", RowBox[{"Exp", "[", RowBox[{ RowBox[{"-", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}]}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], "]"}], " ", RowBox[{ RowBox[{"(", RowBox[{"2", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}], ")"}], "^", "l"}], " ", "*", RowBox[{"LaguerreL", "[", RowBox[{ RowBox[{"n", " ", "-", " ", "l", " ", "-", " ", "1"}], ",", " ", RowBox[{ RowBox[{"2", " ", "l"}], " ", "+", " ", "1"}], ",", " ", RowBox[{"2", " ", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], "*", RowBox[{"Z", "/", RowBox[{"(", RowBox[{"n", "*", "a0"}], ")"}]}]}]}], "]"}], " ", "*", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{ RowBox[{"(", RowBox[{ RowBox[{"2", "l"}], "+", "1"}], ")"}], "/", RowBox[{"(", RowBox[{"2", "Pi"}], ")"}]}], "*", RowBox[{ RowBox[{"Factorial", "[", RowBox[{"l", "-", "m"}], "]"}], "/", RowBox[{"Factorial", "[", RowBox[{"l", "+", "m"}], "]"}]}]}], "]"}], "*", RowBox[{"LegendreP", "[", RowBox[{"l", ",", " ", RowBox[{"Abs", "[", "m", "]"}], ",", " ", RowBox[{"z", "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}]}]}], "]"}], "*", RowBox[{"Cos", "[", RowBox[{"m", "*", RowBox[{"ArcCos", "[", RowBox[{"x", "/", RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}]}], "]"}]}], "]"}]}], "]"}]}], "/;", RowBox[{"m", ">", "0"}]}]}], "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"Psi", "::", "usage"}], "=", "\"\<[{n,l,m},x,y,z] n,l,m are the quantum numbers. x,y,z are the names \ of the x,y and z coordinates\>\""}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", " ", RowBox[{ RowBox[{ RowBox[{"Sqrt", "[", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "]"}], " ", "is", " ", "expressed", " ", "as", " ", RowBox[{"r", ".", " ", "If"}], " ", "you", " ", RowBox[{"don", "'"}], "t", " ", "want", " ", "this", " ", "to", " ", "happen", " ", "delete", " ", "brackets", " ", "after"}], " ", "/.", RowBox[{ RowBox[{"{", " ", "...", " ", "}"}], " ", "in", " ", "the", " ", "following", " ", "expressions"}]}], " ", "*)"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"AtOrbit", " ", "=", " ", RowBox[{"FullSimplify", "[", RowBox[{ RowBox[{ RowBox[{"FullSimplify", "[", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n", ",", "l", ",", "m"}], "}"}], ",", "x", ",", "y", ",", "z"}], "]"}], ",", RowBox[{"Element", "[", RowBox[{ RowBox[{"x", "|", "y", "|", "z"}], ",", "Reals"}], "]"}]}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "\[Rule]", RowBox[{"r", "^", "2"}]}], ",", RowBox[{ RowBox[{"1", "/", RowBox[{"(", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], ")"}]}], "\[Rule]", RowBox[{"1", "/", RowBox[{"r", "^", "2"}]}]}]}], "}"}]}], ",", RowBox[{"r", ">", "0"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{ RowBox[{"LaplaceAtOrbit", " ", "=", " ", RowBox[{"FullSimplify", "[", RowBox[{ RowBox[{ RowBox[{"FullSimplify", "[", RowBox[{ RowBox[{"Laplacian", "[", RowBox[{ RowBox[{"Psi", "[", RowBox[{ RowBox[{"{", RowBox[{"n", ",", "l", ",", "m"}], "}"}], ",", "x", ",", "y", ",", "z"}], "]"}], ",", RowBox[{"{", RowBox[{"x", ",", "y", ",", "z"}], "}"}]}], "]"}], ",", RowBox[{"Element", "[", RowBox[{ RowBox[{"x", "|", "y", "|", "z"}], ",", "Reals"}], "]"}]}], "]"}], "/.", RowBox[{"{", RowBox[{ RowBox[{ RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], "\[Rule]", RowBox[{"r", "^", "2"}]}], ",", RowBox[{ RowBox[{"1", "/", RowBox[{"(", RowBox[{ RowBox[{"x", "^", "2"}], "+", RowBox[{"y", "^", "2"}], "+", RowBox[{"z", "^", "2"}]}], ")"}]}], "\[Rule]", RowBox[{"1", "/", RowBox[{"r", "^", "2"}]}]}]}], "}"}]}], ",", RowBox[{"r", ">", "0"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", "\"\<Atomic Orbital\>\"", "\[IndentingNewLine]", "AtOrbit", "\[IndentingNewLine]", "\[IndentingNewLine]", "\"\<Laplace operator applied to atomic orbital displayed above\>\"", "\[IndentingNewLine]", "LaplaceAtOrbit", "\[IndentingNewLine]", "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"(*", RowBox[{ RowBox[{"Press", " ", "shift"}], " ", "+", RowBox[{"enter", " ", "to", " ", "evaluate"}]}], "*)"}]}]}]], "Input", CellChangeTimes->{ 3.35696210375764*^9, 3.3634725675*^9, {3.680267648969605*^9, 3.6802677777893677`*^9}, {3.680267917025701*^9, 3.680267950118266*^9}, { 3.68026803440273*^9, 3.680268078650187*^9}, {3.680268161660758*^9, 3.680268252117462*^9}, {3.680268283719*^9, 3.680268315847728*^9}, { 3.6802683679885*^9, 3.680268401044072*^9}, {3.6802684347322483`*^9, 3.68026846941199*^9}, 3.68026851589177*^9, {3.680268550847946*^9, 3.680268602485694*^9}, {3.680277200598281*^9, 3.680277200967642*^9}, { 3.680277234870984*^9, 3.680277279999036*^9}, {3.680277368558592*^9, 3.680277379167626*^9}, {3.6802775055031853`*^9, 3.6802775433742323`*^9}, { 3.680278118244079*^9, 3.680278127959669*^9}, {3.6802784585585737`*^9, 3.6802785869469976`*^9}, {3.680278865740156*^9, 3.680278911671702*^9}, { 3.6802791705405188`*^9, 3.680279353827753*^9}, {3.6802794897733517`*^9, 3.6802795064406757`*^9}, {3.680279642796474*^9, 3.6802796473713818`*^9}, { 3.68027969209943*^9, 3.680279693816111*^9}, {3.680279868010215*^9, 3.680279944611951*^9}, {3.680280074808546*^9, 3.680280080471993*^9}, { 3.680280136706254*^9, 3.680280136852335*^9}, {3.680280188871863*^9, 3.680280190484304*^9}, {3.680280326138034*^9, 3.680280332936491*^9}, { 3.68028040865331*^9, 3.680280409902307*^9}, {3.680422742828141*^9, 3.680422743874589*^9}, {3.680454127131035*^9, 3.680454139436511*^9}, { 3.6804542844683437`*^9, 3.680454507163357*^9}, {3.6804982821582613`*^9, 3.6804983117917747`*^9}, {3.680498545809291*^9, 3.6804986280861063`*^9}, { 3.680498859127798*^9, 3.6804990405705233`*^9}, {3.680499458218275*^9, 3.6804995193350782`*^9}, {3.680499994141029*^9, 3.680500114217901*^9}, { 3.680500145428303*^9, 3.6805001801054068`*^9}, {3.680500356902849*^9, 3.680500381139038*^9}, {3.6805004478266573`*^9, 3.68050070843865*^9}, { 3.68050080824401*^9, 3.680500830439155*^9}, {3.6805008797742367`*^9, 3.68050096876678*^9}, {3.6805010324238863`*^9, 3.680501038408432*^9}, { 3.680501083067733*^9, 3.6805011033392887`*^9}, {3.680501473869156*^9, 3.680501475939424*^9}, {3.680501507238729*^9, 3.680501534137072*^9}, { 3.680501760785087*^9, 3.680501840435904*^9}, {3.680501878359671*^9, 3.6805018960181627`*^9}, {3.680536350151812*^9, 3.680536420825638*^9}, 3.680536465213065*^9, {3.680536659766471*^9, 3.6805366985140343`*^9}, { 3.680536737643167*^9, 3.680536743236704*^9}, {3.680536829922267*^9, 3.6805368304467*^9}, {3.6805369171100407`*^9, 3.6805370292866*^9}, { 3.680537116888463*^9, 3.6805371337379723`*^9}, 3.6805371761775227`*^9}], Cell[BoxData["\<\"Atomic Orbital\"\>"], "Output", CellChangeTimes->{{3.680499498899252*^9, 3.680499520058292*^9}, 3.6805001154586563`*^9, {3.680500155454672*^9, 3.680500181483968*^9}, 3.6805006880295258`*^9, {3.680500943444726*^9, 3.6805009694356823`*^9}, 3.6805010395271378`*^9, 3.680501087404969*^9, 3.680501476918788*^9, { 3.680501509067039*^9, 3.680501534991123*^9}, 3.680501781348852*^9, 3.680501818754817*^9, 3.680501849812171*^9, {3.680536400478294*^9, 3.6805364238221607`*^9}, {3.680536688791065*^9, 3.680536700710106*^9}, 3.680536925115191*^9, {3.680536956377844*^9, 3.680536978694578*^9}, 3.680537032284369*^9, 3.680537412786456*^9}], Cell[BoxData[ FractionBox[ RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"-", FractionBox[ RowBox[{"r", " ", "Z"}], RowBox[{"2", " ", "a0"}]]}]], " ", "z", " ", "Z", " ", SqrtBox[ FractionBox[ SuperscriptBox["Z", "3"], SuperscriptBox["a0", "3"]]]}], RowBox[{"4", " ", "a0", " ", SqrtBox[ RowBox[{"2", " ", "\[Pi]"}]]}]]], "Output", CellChangeTimes->{{3.680499498899252*^9, 3.680499520058292*^9}, 3.6805001154586563`*^9, {3.680500155454672*^9, 3.680500181483968*^9}, 3.6805006880295258`*^9, {3.680500943444726*^9, 3.6805009694356823`*^9}, 3.6805010395271378`*^9, 3.680501087404969*^9, 3.680501476918788*^9, { 3.680501509067039*^9, 3.680501534991123*^9}, 3.680501781348852*^9, 3.680501818754817*^9, 3.680501849812171*^9, {3.680536400478294*^9, 3.6805364238221607`*^9}, {3.680536688791065*^9, 3.680536700710106*^9}, 3.680536925115191*^9, {3.680536956377844*^9, 3.680536978694578*^9}, 3.680537032284369*^9, 3.680537412787899*^9}], Cell[BoxData["\<\"Laplace operator applied to atomic orbital displayed \ above\"\>"], "Output", CellChangeTimes->{{3.680499498899252*^9, 3.680499520058292*^9}, 3.6805001154586563`*^9, {3.680500155454672*^9, 3.680500181483968*^9}, 3.6805006880295258`*^9, {3.680500943444726*^9, 3.6805009694356823`*^9}, 3.6805010395271378`*^9, 3.680501087404969*^9, 3.680501476918788*^9, { 3.680501509067039*^9, 3.680501534991123*^9}, 3.680501781348852*^9, 3.680501818754817*^9, 3.680501849812171*^9, {3.680536400478294*^9, 3.6805364238221607`*^9}, {3.680536688791065*^9, 3.680536700710106*^9}, 3.680536925115191*^9, {3.680536956377844*^9, 3.680536978694578*^9}, 3.680537032284369*^9, 3.680537412788458*^9}], Cell[BoxData[ FractionBox[ RowBox[{ SuperscriptBox["\[ExponentialE]", RowBox[{"-", FractionBox[ RowBox[{"r", " ", "Z"}], RowBox[{"2", " ", "a0"}]]}]], " ", "z", " ", SuperscriptBox["Z", "2"], " ", SqrtBox[ FractionBox[ SuperscriptBox["Z", "3"], RowBox[{ SuperscriptBox["a0", "3"], " ", SuperscriptBox["r", "2"]}]]], " ", RowBox[{"(", RowBox[{ RowBox[{ RowBox[{"-", "8"}], " ", "a0"}], "+", RowBox[{"r", " ", "Z"}]}], ")"}]}], RowBox[{"16", " ", SuperscriptBox["a0", "3"], " ", SqrtBox[ RowBox[{"2", " ", "\[Pi]"}]]}]]], "Output", CellChangeTimes->{{3.680499498899252*^9, 3.680499520058292*^9}, 3.6805001154586563`*^9, {3.680500155454672*^9, 3.680500181483968*^9}, 3.6805006880295258`*^9, {3.680500943444726*^9, 3.6805009694356823`*^9}, 3.6805010395271378`*^9, 3.680501087404969*^9, 3.680501476918788*^9, { 3.680501509067039*^9, 3.680501534991123*^9}, 3.680501781348852*^9, 3.680501818754817*^9, 3.680501849812171*^9, {3.680536400478294*^9, 3.6805364238221607`*^9}, {3.680536688791065*^9, 3.680536700710106*^9}, 3.680536925115191*^9, {3.680536956377844*^9, 3.680536978694578*^9}, 3.680537032284369*^9, 3.680537412789257*^9}] }, Open ]], Cell[CellGroupData[{ Cell[TextData[{ StyleBox["Evaluate following cells to get ", "Chapter"], Cell[BoxData[ FormBox[ StyleBox[ RowBox[{"L", StyleBox[ AdjustmentBox["A", BoxBaselineShift->-0.4, BoxMargins->{{-0.5, -0.3}, {0, 0}}], FontSize->Smaller], "T", AdjustmentBox["E", BoxBaselineShift->0.5, BoxMargins->{{-0.3, 0}, {0, 0}}], "X"}], SingleLetterItalics->False], TraditionalForm]]], StyleBox[" syntax", "Chapter"], " " }], "Section", CellChangeTimes->{{3.680537047447669*^9, 3.6805370920176907`*^9}, { 3.68053741912608*^9, 3.680537449967037*^9}, {3.680537552860992*^9, 3.680537560032627*^9}}], Cell[CellGroupData[{ Cell[BoxData[{"\"\<Atomic Orbital\>\"", "\[IndentingNewLine]", RowBox[{"TeXForm", "[", "AtOrbit", "]"}]}], "Input", CellChangeTimes->{{3.680537452108548*^9, 3.680537509094211*^9}, { 3.680537548839533*^9, 3.6805375672542763`*^9}}], Cell[BoxData["\<\"Atomic Orbital\"\>"], "Output", CellChangeTimes->{3.6805375694238787`*^9}], Cell["\<\ \\frac{z Z \\sqrt{\\frac{Z^3}{\\text{a0}^3}} e^{-\\frac{r Z}{2 \\text{a0}}}}{4 \\sqrt{2 \\pi } \\text{a0}}\ \>", "Output", CellChangeTimes->{3.68053756942557*^9}] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{"\[IndentingNewLine]", RowBox[{ "\"\<Laplace operator applied to atomic orbital displayed above\>\"", "\[IndentingNewLine]", RowBox[{"TeXForm", "[", "LaplaceAtOrbit", "]"}]}]}]], "Input", CellChangeTimes->{{3.680537452108548*^9, 3.680537509094211*^9}, { 3.680537548839533*^9, 3.6805375672542763`*^9}}], Cell[BoxData["\<\"Laplace operator applied to atomic orbital displayed \ above\"\>"], "Output", CellChangeTimes->{3.680537511655724*^9, 3.6805375723269787`*^9}], Cell["\<\ \\frac{z Z^2 e^{-\\frac{r Z}{2 \\text{a0}}} \\sqrt{\\frac{Z^3}{\\text{a0}^3 r^2}} (r Z-8 \\text{a0})}{16 \\sqrt{2 \\pi } \\text{a0}^3}\ \>", "Output", CellChangeTimes->{3.680537511655724*^9, 3.680537572328401*^9}] }, Open ]] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["\<\ Evaluate following cells to get orbitals in mathematica syntax\ \>", "Chapter", CellChangeTimes->{{3.6805375953984947`*^9, 3.680537709515551*^9}}], Cell[CellGroupData[{ Cell[BoxData[{"\"\<Atomic Orbital\>\"", "\[IndentingNewLine]", RowBox[{"InputForm", "[", "AtOrbit", "]"}]}], "Input", CellChangeTimes->{{3.6805375953984947`*^9, 3.6805376593171663`*^9}, 3.6805377145532017`*^9}], Cell[BoxData["\<\"Atomic Orbital\"\>"], "Output", CellChangeTimes->{3.6805377176185207`*^9}], Cell["(z*Z*Sqrt[Z^3/a0^3])/(4*a0*E^((r*Z)/(2*a0))*Sqrt[2*Pi])", "Output", CellChangeTimes->{3.6805377176225777`*^9}] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[{"\"\<Laplace operator applied to atomic orbital displayed above\ \>\"", "\[IndentingNewLine]", RowBox[{"InputForm", "[", "LaplaceAtOrbit", "]"}]}], "Input", CellChangeTimes->{{3.6805375953984947`*^9, 3.6805376593171663`*^9}, 3.6805377145532017`*^9}], Cell[BoxData["\<\"Laplace operator applied to atomic orbital displayed \ above\"\>"], "Output", CellChangeTimes->{3.680537719968749*^9}], Cell["\<\ (z*Z^2*Sqrt[Z^3/(a0^3*r^2)]*(-8*a0 + r*Z))/ (16*a0^3*E^((r*Z)/(2*a0))*Sqrt[2*Pi])\ \>", "Output", CellChangeTimes->{3.6805377199694843`*^9}] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell["Calculate overlap integral", "Chapter", CellChangeTimes->{{3.680501456845961*^9, 3.680501462062808*^9}}], Cell[CellGroupData[{ Cell[BoxData[ RowBox[{ RowBox[{"(*", " ", RowBox[{ RowBox[{ "If", " ", "you", " ", "want", " ", "to", " ", "calculate", " ", "an", " ", "overlap", " ", "integral"}], ",", " ", RowBox[{ "evaluate", " ", "following", " ", "cell", " ", "and", " ", "specify", " ", "quantum", " ", "numbers", " ", "and", " ", "distance", " ", "between", " ", "atoms"}]}], "*)"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{ RowBox[{ RowBox[{"a0", " ", "=", " ", RowBox[{"5.291772106712", "*", RowBox[{"10", "^", RowBox[{"-", "11"}]}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"Z", "=", "1"}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"d", " ", "=", " ", RowBox[{"74", "*", RowBox[{"10", "^", RowBox[{"-", "12"}]}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"quantumnumbers1", "=", RowBox[{"{", RowBox[{"1", ",", "0", ",", "0"}], "}"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"quantumnumbers2", "=", RowBox[{"{", RowBox[{"1", ",", "0", ",", "0"}], "}"}]}], ";"}], "\[IndentingNewLine]", "\[IndentingNewLine]", RowBox[{"NIntegrate", "[", RowBox[{ RowBox[{ RowBox[{"Psi", "[", RowBox[{"quantumnumbers1", ",", RowBox[{"x", "+", RowBox[{"d", "/", "2"}]}], ",", "y", ",", "z"}], "]"}], "*", RowBox[{"Psi", "[", RowBox[{"quantumnumbers2", ",", RowBox[{"x", "-", RowBox[{"d", "/", "2"}]}], ",", "y", ",", "z"}], "]"}]}], ",", RowBox[{"{", RowBox[{"x", ",", RowBox[{ RowBox[{"-", "50"}], "a0"}], ",", RowBox[{"50", "a0"}]}], "}"}], ",", RowBox[{"{", RowBox[{"y", ",", RowBox[{ RowBox[{"-", "50"}], "a0"}], ",", RowBox[{"50", "a0"}]}], "}"}], ",", RowBox[{"{", RowBox[{"z", ",", RowBox[{ RowBox[{"-", "50"}], "a0"}], ",", RowBox[{"50", "a0"}]}], "}"}], ",", RowBox[{"WorkingPrecision", "\[Rule]", "50"}]}], "]"}]}]}]], "Input", CellChangeTimes->{{3.680419309320182*^9, 3.6804194126041203`*^9}, { 3.680419526994852*^9, 3.68041960922066*^9}, {3.680419988741125*^9, 3.680420004703898*^9}, {3.6804203907187023`*^9, 3.68042041703402*^9}, { 3.6804204879400177`*^9, 3.68042050686876*^9}, {3.680420577425701*^9, 3.6804206821812553`*^9}, {3.6804208435592623`*^9, 3.680420878992793*^9}, { 3.680421166685328*^9, 3.6804212211179237`*^9}, 3.680421477851687*^9, { 3.680421543597632*^9, 3.680421617862692*^9}, {3.680421653005362*^9, 3.680421656577001*^9}, {3.680421744649856*^9, 3.680421745727749*^9}, { 3.680421975341563*^9, 3.680421994260213*^9}, {3.680422137026458*^9, 3.680422142286463*^9}, {3.680422572513933*^9, 3.680422573066861*^9}, { 3.680422677741963*^9, 3.68042271435723*^9}, {3.68042278698472*^9, 3.6804227873066807`*^9}, {3.6804229545240088`*^9, 3.6804229913791113`*^9}, {3.680437926126237*^9, 3.680437929634451*^9}, { 3.68050026284394*^9, 3.680500266092196*^9}, {3.680501167413177*^9, 3.6805013432373657`*^9}}], Cell[BoxData["0.\ 75338532147689655474447162064044959127809752867891332473544597952866095782077`\ 50."], "Output", CellChangeTimes->{3.680501408372034*^9}] }, Open ]] }, Open ]] }, WindowSize->{1280, 1376}, WindowMargins->{{-1, Automatic}, {Automatic, -1}}, Magnification->1.5, FrontEndVersion->"10.4 for Linux x86 (64-bit) (April 11, 2016)", StyleDefinitions->"Default.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[558, 20, 2402, 40, 234, "Input"], Cell[CellGroupData[{ Cell[2985, 64, 1999, 27, 98, "Chapter"], Cell[CellGroupData[{ Cell[5009, 95, 16187, 415, 1445, "Input"], Cell[21199, 512, 675, 9, 45, "Output"], Cell[21877, 523, 1022, 23, 131, "Output"], Cell[22902, 548, 721, 10, 45, "Output"], Cell[23626, 560, 1262, 32, 131, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[24925, 597, 648, 20, 109, "Section"], Cell[CellGroupData[{ Cell[25598, 621, 234, 3, 75, "Input"], Cell[25835, 626, 93, 1, 45, "Output"], Cell[25931, 629, 176, 4, 87, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[26144, 638, 340, 7, 104, "Input"], Cell[26487, 647, 161, 2, 45, "Output"], Cell[26651, 651, 230, 5, 111, "Output"] }, Open ]] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[26942, 663, 157, 3, 165, "Chapter"], Cell[CellGroupData[{ Cell[27124, 670, 217, 3, 75, "Input"], Cell[27344, 675, 93, 1, 45, "Output"], Cell[27440, 678, 117, 1, 62, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[27594, 684, 270, 4, 75, "Input"], Cell[27867, 690, 137, 2, 45, "Output"], Cell[28007, 694, 152, 4, 87, "Output"] }, Open ]] }, Open ]], Cell[CellGroupData[{ Cell[28208, 704, 111, 1, 98, "Chapter"], Cell[CellGroupData[{ Cell[28344, 709, 3135, 76, 336, "Input"], Cell[31482, 787, 156, 3, 45, "Output"] }, Open ]] }, Open ]] } ] *) function [ result ] = atorbit(xlocal,ylocal,zlocal,Zlocal,laplace,type) %ATORBIT Atomic orbital centered at specified coordinates % % xlocal ... x-coordinate of the center of the atom (in m) % ylocal ... y-coordinate (in m) % zlocal ... z-coordinate (in m) % Zlocal ... atomic number of the atom % laplace ... boolean. true - laplace operator is appled to atomic orbital % false - "pure" atomic orbital is returned % type ... string specifying the orbital (1s, 2s, 2px, 2py, 2pz) result = []; %[x_mesh,xAtom_mesh] = meshgrid(x,xAtom); %xlocal = x - xAtom; %[y_mesh,yAtom_mesh] = meshgrid(y,yAtom); %ylocal = y - yAtom; %[z_mesh,zAtom_mesh] = meshgrid(z,zAtom); %zlocal = z - zAtom; %Zlocal = repmat(Z,length(x),1); a0 = 5.2917721E-11; % Bohr radius in meters r = sqrt(xlocal.^2 + ylocal.^2 + zlocal.^2); if ~laplace switch lower(type) case '1s' result = sqrt(Zlocal.^3./(pi*a0.^3)).*exp(-Zlocal.*r/a0); case '2s' result= 1/4*sqrt(Zlocal.^3./(2*pi*a0.^3)).*(2-Zlocal.*r/a0).*exp(-Zlocal.*r/(2*a0)); case '2px' result= 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*xlocal.*exp(-Zlocal.*r/(2*a0)); case '2py' result= 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*ylocal.*exp(-Zlocal.*r/(2*a0)); case '2pz' result= 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*zlocal.*exp(-Zlocal.*r/(2*a0)); end else switch(lower(type)) case '1s' result = sqrt(Zlocal.^5./(pi*a0.^5)).*(Zlocal/a0 - 2./r).*exp(-Zlocal.*r/a0); case '2s' result = 1/4*sqrt(Zlocal.^5./(2*pi*a0.^5)).*(-Zlocal.^2.*r/(4*a0^2) + 5*Zlocal/(2*a0) - 4./r).*exp(-Zlocal.*r/(2*a0)); case '2px' result = 1/4*sqrt(Zlocal.^7./(2*pi*a0.^7)).*(Zlocal.*r/(4*a0)- 2).*xlocal./r.*exp(-Zlocal.*r/(2*a0)); case '2py' result = 1/4*sqrt(Zlocal.^7./(2*pi*a0.^7)).*(Zlocal.*r/(4*a0)- 2).*ylocal./r.*exp(-Zlocal.*r/(2*a0)); case '2pz' result = 1/4*sqrt(Zlocal.^7./(2*pi*a0.^7)).*(Zlocal.*r/(4*a0)- 2).*zlocal./r.*exp(-Zlocal.*r/(2*a0)); end end end #include <stdio.h> #include <stdlib.h> #include <conio.h> #include <math.h> #include <time.h> #define A0 ( 1 ) #define HA_EV ( 27.211386 ) /** Pointer to atomic orbital function x. */ typedef double (*ao_t)(double x, double y, double z, double xi, double yi, double zi, double Za); /** Pointer to Laplacian of atomic orbital x. */ typedef double (*L_ao_t)(double x, double y, double z, double xi, double yi, double zi, double Za); double ao_1s(double x, double y, double z, double xi, double yi, double zi, double Za); double ao_2s(double x, double y, double z, double xi, double yi, double zi, double Za); double ao_2px(double x, double y, double z, double xi, double yi, double zi, double Za); double ao_2py(double x, double y, double z, double xi, double yi, double zi, double Za); double ao_2pz(double x, double y, double z, double xi, double yi, double zi, double Za); double ao_3s(double x, double y, double z, double xi, double yi, double zi, double Za); double L_ao_1s(double x, double y, double z, double xi, double yi, double zi, double Za); double L_ao_2s(double x, double y, double z, double xi, double yi, double zi, double Za); double L_ao_2px(double x, double y, double z, double xi, double yi, double zi, double Za); double L_ao_2py(double x, double y, double z, double xi, double yi, double zi, double Za); double L_ao_2pz(double x, double y, double z, double xi, double yi, double zi, double Za); double L_ao_3s(double x, double y, double z, double xi, double yi, double zi, double Za); /** Returns a random number uniformly distributed between 0 and 1. */ static double getRandom(void) { return ((double)rand()/(double)RAND_MAX); } /** Calculates the radius from given Cartesian coordinates. */ static double getRadius(double x, double y, double z) { return sqrt(pow(x,2)+pow(y,2)+pow(z,2)); } int main() { int i; double Za, E, x_a, x_b, y_a, y_b, z_a, z_b, xi_a, xi_b, yi_a, yi_b, zi_a, zi_b; ao_t ao; L_ao_t L_ao; // feed the random number generator with the system time as seed srand((unsigned)time(NULL)); /* Atomic units are used at this point to simplify the Hamiltonian: H = - nabla^2/2 - Za/r Physical quantities are therefore given with units different from SI: - length: in terms of Bohr-lengths, 1 A0 = 0.52918 Angstroem - energy: unit Hartree, 1 Ha = 2 Ry = 27.211 eV To calculate the energy we start from the equation H phi = E phi, and convert it to E = (H phi) / phi; with the atomic orbitals phi = ao and the related Laplacian ( nabla^2 ao ) = L_ao we obtain E = -0.5*L_ao/ao - Za/r; This will be done for different random positions (xi-x,yi-y,zi-z) for the electron of the atom (xi,yi,zi) to confirm the atomic orbital function and the related Laplacian function. If the energy stays the same for different positions the given atomic orbital can be assumed to be confirmed. */ /* ********** select parameters ***************************************** */ // atomic number Za = 1.0; // range of random numbers x_a = -20; x_b = 20; y_a = -20; y_b = 20; z_a = -10; z_b = 10; xi_a = -50; xi_b = 50; yi_a = -50; yi_b = 50; zi_a = 0; zi_b = 100; // select an atomic orbital function to be tested ao = &ao_3s; L_ao = &L_ao_3s; /* ********************************************************************** */ printf("atomic number: %d\n\n", (int)Za); // calculate energy for 10 random positions for(i=0; i<10; i++) { double x, y, z, xi, yi, zi, r; x = getRandom()*(x_b-x_a)+x_a; y = getRandom()*(y_b-y_a)+y_a; z = getRandom()*(z_b-z_a)+z_a; xi = getRandom()*(xi_b-xi_a)+xi_a; yi = getRandom()*(yi_b-yi_a)+yi_a; zi = getRandom()*(zi_b-zi_a)+zi_a; r = getRadius(x-xi, y-yi, z-zi); E = -0.5*L_ao(x,y,z,xi,yi,zi,Za)/ao(x,y,z,xi,yi,zi,Za) - Za/r; // note: double has about 16 significant figures (53 bits of precision) printf("E: %8.6lf Ha = %8.6lf eV\n", E, E*HA_EV); } // wait for any key pressed getch(); return 0; } double ao_1s(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return sqrt(pow(Za,3)/(M_PI*pow(A0,3))) * exp(-Za*r/A0); } double ao_2s(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return 0.25 * sqrt(pow(Za,3)/(2.*M_PI*pow(A0,3))) * (2. - Za*r/A0) * exp(-Za*r/(2.*A0)); } double ao_2px(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return 0.25 * sqrt(pow(Za,5)/(2.*M_PI*pow(A0,5))) * (x - xi) * exp(-Za*r/(2.*A0)); } double ao_2py(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return 0.25 * sqrt(pow(Za,5)/(2.*M_PI*pow(A0,5))) * (y - yi) * exp(-Za*r/(2.*A0)); } double ao_2pz(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return 0.25 * sqrt(pow(Za,5)/(2.*M_PI*pow(A0,5))) * (z - zi) * exp(-Za*r/(2.*A0)); } double ao_3s(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return 1./81 * sqrt(pow(Za,3)/(3.*M_PI*pow(A0,3))) * (27. - 18.*Za*r/A0 + 2.*pow(Za,2)*pow(r,2)/pow(A0,2)) * exp(-Za*r/(3.*A0)); } double L_ao_1s(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return sqrt(pow(Za,5)/(M_PI*pow(A0,5))) * (Za/A0 - 2./r) * exp(-Za*r/A0); } double L_ao_2s(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return 0.25 * sqrt(pow(Za,5)/(2.*M_PI*pow(A0,5))) * (-Za*Za*r/(4.*pow(A0,2)) + 5.*Za/(2*A0) - 4./r) * exp(-Za*r/(2.*A0)); } double L_ao_2px(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return sqrt(pow(Za,7)/(pow(2,5)*M_PI*pow(A0,7))) * (Za*r/(4.*A0) - 2.) * ((x - xi)/r) * exp(-Za*r/(2.*A0)); } double L_ao_2py(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return sqrt(pow(Za,7)/(pow(2,5)*M_PI*pow(A0,7))) * (Za*r/(4.*A0) - 2.) * ((y - yi)/r) * exp(-Za*r/(2.*A0)); } double L_ao_2pz(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return sqrt(pow(Za,7)/(pow(2,5)*M_PI*pow(A0,7))) * (Za*r/(4.*A0) - 2.) * ((z - zi)/r) * exp(-Za*r/(2.*A0)); } double L_ao_3s(double x, double y, double z, double xi, double yi, double zi, double Za) { double r = getRadius(x-xi, y-yi, z-zi); return sqrt(pow(Za,5)/(3.*M_PI*pow(A0,5))) * (-486. + 351.*r*Za/A0 - 54.*pow(Za,2)*pow(r,2)/(pow(A0,2)) + 2.*pow(Za,3)*pow(r,3)/pow(A0,3)) * exp(-Za*r/(3.*A0)) / (729.*r); } // EOF <script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = 3.1415926535897932; var hbar = 1.05457266E-34; var ec = 1.60217733E-19; var eps0 = 8.854187817E-12; var me = 9.1093897E-31; var k1 = 0.5*hbar*hbar/me; var k2 = ec*ec/(4*pi*eps0); //Coulomb constant function ao_1s(x,y,z,xi,yi,zi,Za) { //the 1s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.sqrt(Math.pow(Za,3)/(pi*Math.pow(a0,3)))*Math.exp(-Za*r/a0); } function L_ao_1s(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 1s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.sqrt(Math.pow(Za,5)/(pi*Math.pow(a0,5)))*(Za/a0 -2/r)*Math.exp(-Za*r/a0); } function ao_2s(x,y,z,xi,yi,zi,Za) { //the 2s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return 0.25*Math.sqrt(Math.pow(Za,3)/(2*pi*Math.pow(a0,3)))*(2-Za*r/a0)*Math.exp(-Za*r/(2*a0)); } function L_ao_2s(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 2s atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return 0.25*Math.sqrt(Math.pow(Za,5)/(2*pi*Math.pow(a0,5)))*(-Za*Za*r/(4*a0*a0)+5*Za/(2*a0)-4/r)*Math.exp(-Za*r/(2*a0)); } function ao_2pz(x,y,z,xi,yi,zi,Za) { //the 2pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.sqrt(Math.pow(Za/(2*a0),5)/pi)*(z-zi)*Math.exp(-Za*r/(2*a0)); } function L_ao_2pz(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 2pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2)+Math.pow(y-yi,2)+Math.pow(z-zi,2)); return Math.pow(Za/a0,7/2)*Math.pow(1/2,5/2)*Math.sqrt(1/pi)*(Za*r/(4*a0)-2)*((z-zi)/r)*Math.exp(-Za*r/(2*a0)); } function ao_3s(x,y,z,xi,yi,zi,Za) { //the 3s atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(Math.pow(Za/a0,3) / (3*pi)) * (27 - 18*Za/a0*r + 2*Math.pow(Za/a0*r, 2)) * Math.exp(-Za/a0 * r/3); } function L_ao_3s(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3s atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(Math.pow(Za/a0,3) / (3*pi)) * (-54*Za/a0/r + 39*Math.pow(Za/a0,2) - 6*Math.pow(Za/a0,3)*r + 2/9*Math.pow(Za/a0,4)*r*r) * Math.exp(-Za/a0 * r/3); } function ao_3px(x,y,z,xi,yi,zi,Za) { //the 3px atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (6 - Za/a0*r) * Za/a0*(x-xi) * Math.exp(-Za/a0 * r/3); } function L_ao_3px(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3px atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (-12*Za/a0/r + 8/3*Math.pow(Za/a0,2) - Math.pow(Za/a0,3)*r/9) * Za/a0*(x-xi) * Math.exp(-Za/a0 * r/3); } function ao_3py(x,y,z,xi,yi,zi,Za) { //the 3py atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (6 - Za/a0*r) * Za/a0*(y-yi) * Math.exp(-Za/a0 * r/3); } function L_ao_3py(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3py atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (-12*Za/a0/r + 8/3*Math.pow(Za/a0,2) - Math.pow(Za/a0,3)*r/9) * Za/a0*(y-yi) * Math.exp(-Za/a0 * r/3); } function ao_3pz(x,y,z,xi,yi,zi,Za) { //the 3pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (6 - Za/a0*r) * Za/a0*(z-zi) * Math.exp(-Za/a0 * r/3); } function L_ao_3pz(x,y,z,xi,yi,zi,Za) { //Lagrangian of the 3pz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2 * Math.pow(Za/a0,3) / pi) * (-12*Za/a0/r + 8/3*Math.pow(Za/a0,2) - Math.pow(Za/a0,3)*r/9) * Za/a0*(z-zi) * Math.exp(-Za/a0 * r/3); } function ao_3dxz(x,y,z,xi,yi,zi,Za) { //the 3dxz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2*Za**7/(pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((x-xi)*(z-zi))/r**2; } function L_ao_3dxz(x,y,z,xi,yi,zi,Za) { //the 3dxz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(2*Za**9/(pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((x-xi)*(z-zi))/r**2; } function ao_3dxy(x,y,z,xi,yi,zi,Za) { //the 3dxy atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2*Za**7/(pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((x-xi)*(y-yi))/r**2; } function L_ao_3dxy(x,y,z,xi,yi,zi,Za) { //the 3dxy atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(2*Za**9/(pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((x-xi)*(y-yi))/r**2; } function ao_3dyz(x,y,z,xi,yi,zi,Za) { //the 3dyz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(2*Za**7/(pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((y-yi)*(z-zi))/r**2; } function L_ao_3dyz(x,y,z,xi,yi,zi,Za) { //the 3dyz atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(2*Za**9/(pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((y-yi)*(z-zi))/r**2; } function ao_3dzz(x,y,z,xi,yi,zi,Za) { //the 3dz^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81*Math.sqrt(Za**7/(6*pi*a0**7))*Math.exp(-Za*r/(3*a0))*r**2*(3*(z-zi)**2/r**2-1); } function L_ao_3dzz(x,y,z,xi,yi,zi,Za) { //the 3dz^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243*Math.sqrt(Za**9/(6*pi*a0**9))*Math.exp(-Za*r/(3*a0))*r*(Za*r/(3*a0)-6)*(3*(z-zi)**2/r**2-1); } function ao_3dx2my2(x,y,z,xi,yi,zi,Za) { //the 3dx^2-y^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/81 * Math.sqrt(Za**7/(2*pi*a0**7))*r**2*Math.exp(-Za*r/(3*a0))*((x-xi)*(x-xi)-(y-yi)*(y-yi))/r**2; } function L_ao_3dx2my2(x,y,z,xi,yi,zi,Za) { //the 3dx^2-y^2 atomic orbital r = Math.sqrt(Math.pow(x-xi,2) + Math.pow(y-yi,2) + Math.pow(z-zi,2)); return 1/243 * Math.sqrt(Za**9/(2*pi*a0**9))*r*(Za*r/(3*a0)-6)*Math.exp(-Za*r/(3*a0))*((x-xi)*(x-xi)-(y-yi)*(y-yi))/r**2; } aos = [ao_1s, ao_2s, ao_2pz, ao_3s, ao_3px, ao_3py, ao_3pz, ao_3dzz, ao_3dxy, ao_3dxz, ao_3dyz, ao_3dx2my2]; L_aos = [L_ao_1s, L_ao_2s, L_ao_2pz, L_ao_3s, L_ao_3px, L_ao_3py, L_ao_3pz, L_ao_3dzz, L_ao_3dxy, L_ao_3dxz, L_ao_3dyz, L_ao_3dx2my2]; Za = 1; for(j=0; j<aos.length; j++) { ao = aos[j]; L_ao = L_aos[j]; // Print orbital name document.write(ao.name + "<br>"); // Integrate density to check normalization sum = 0; N = 100000; A = 20 * a0; // Half side length of cube over which should be integrated for (i=0; i<N; i++) { // Monte-Carlo integration x = (2*Math.random()-1) * A; xi = (2*Math.random()-1) * A/4; y = (2*Math.random()-1) * A; yi = (2*Math.random()-1) * A/4; z = (2*Math.random()-1) * A; zi = (2*Math.random()-1) * A/4; sum += Math.pow(ao(x, y, z, xi, yi, zi, Za), 2); } document.write("Normalization: "+Math.pow(2*A, 3) * sum/N + "<br>"); // Calculate energy at random points for (i=0; i<10; i++) { x = (2*Math.random()-1) * a0; xi = (2*Math.random()-1) * a0; y = (2*Math.random()-1) * a0; yi = (2*Math.random()-1) * a0; z = (2*Math.random()-1) * a0; zi = (2*Math.random()-1) * a0; E = -k1*L_ao(x,y,z,xi,yi,zi,Za)/ao(x,y,z,xi,yi,zi,Za)-k2*Za/r; document.write(E/ec+" eV<br>"); } document.write("<br>"); } </script>

The following code checks the energies and the normalizations of the atomic orbitals up to 3d.

---

• A Jupyter notebook for plotting the atomic orbitals