1s Atomic Orbital

The 1s atomic orbital is,

\[\begin{equation} \phi^Z_{1s}=\sqrt{\frac{Z^3}{\pi a^3_0}}\exp\left(-\frac{Zr}{a_0}\right). \end{equation}\]

$a_0^{3/2}\phi_{1s}$
	$\frac{x}{a_0}$

To check that the 1s orbital solves the Schrödinger 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}

we let the Laplacian operator act on the orbital,

\[\begin{equation} \nabla^2\phi^Z_{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), \end{equation}\]

and then calculate the energy $E_1= - \frac{\hbar^2}{2m}\frac{\nabla^2\phi_{1s}(\vec{r})}{\phi_{1s}(\vec{r})} - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}|}$. The code below defines the 1s orbital and its Laplacian in Cartesian coordinates centered at position $(x_i,y_i,z_i)$. It then chooses random numbers for $x$, $y$, $z$, $x_i$, $y_i$, and $y_i$ and calculates the energy. If the orbital has been programmed properly, the energy should be -13.6 eV for every position.

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); 
}

for (i=0; i<10; i++) { //The energy is calculated for 10 random positions
  x = 2*a0*(0.5-Math.random())/Z_eff;
  xi = 2*a0*(0.5-Math.random())/Z_eff;
  y = 2*a0*(0.5-Math.random())/Z_eff;
  yi = 2*a0*(0.5-Math.random())/Z_eff;
  z = 2*a0*(0.5-Math.random())/Z_eff;
  zi = 2*a0*(0.5-Math.random())/Z_eff;
  E = -k1*L_ao_1s(x,y,z,xi,yi,zi,Z_eff)/ao_1s(x,y,z,xi,yi,zi,Z_eff)-k2*Z_eff/r;
  document.write("E<sub>1</sub> = "+E/ec+" eV<br>");
}
</script>

The probability of finding a electron a distance $r$ from the nucleus is $P(r)=4\pi r^2|\phi_{1s}|^2$. The probability has a maximum at $a_0$ but by looking at the integral is clear that it is more probable to find the electrons further than $a_0$ from the nucleus than closer than $a_0$ from the nucleus.


	$\frac{Z_{eff}r}{a_0}$

Consider the overlap integral of two 1s orbitals located at positions $\vec{r}_1$ and $\vec{r}_2$,

\begin{equation} S_{12}=\langle \phi_{1s}(\vec{r}-\vec{r}_1)|\phi_{1s}(\vec{r}-\vec{r}_2)\rangle . \end{equation}

For a hydrogen molecule, $\vec{r}_1=-0.38\,\hat{x}$ Å and $\vec{r}_2=0.38\,\hat{x}$ Å. Below $\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)$ is plotted along the $x$-axis and along the $y$-axis.

$a_0^3\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)$
	$\frac{x}{a_0}\,,\frac{y}{a_0}$

The code below uses a Monte-Carlo method to integrate $\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)$ and calculate $S_{12}$.

<script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = Math.PI; var Z_eff = 1; 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)/(piMath.pow(a0,3)))Math.exp(-Zar/a0); } S12 = 0; N = 1000000; //number of random points x1 = -0.38E-10; x2 = 0.38E-10; L = 8a0/Z_eff for (i=0; i<N; i++) { x = L(0.5-Math.random()); y = L(0.5-Math.random()); z = L(0.5-Math.random()); S12 += ao_1s(x,y,z,x1,0,0,Z_eff)ao_1s(x,y,z,x2,0,0,Z_eff); } document.write("S<sub>12<\/sub> = "+Math.pow(L,3)*S12/N); //multiply by the volume and divide by the number of points </script>

Press the 'Execute' button a few times and notice that the answer keeps changing. By doing this you can estimate the error in the calculation. The error should decrease like $1/\sqrt{N}$ where $N$ is the number of random numbers chosen. By setting $x_1=x_2=0$ in the code and pressing 'Execute', you calculate $\langle \phi_{1s}(\vec{r})|\phi_{1s}(\vec{r})\rangle$ which should equal 1 if the wave function is properly normalized.

The diagonal Hamiltonian matrix element of a homonuclear diatomic molecule (H₂, O₂, N₂, etc.) with two 1s orbitals located at position $\vec{r}_1$ is,

\begin{equation} H_{11}=\Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}- \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}^Z(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

This can be broken into two terms,

\begin{equation} H_{11}=\Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}\right|\phi_{1s}^Z(\vec{r}-\vec{r}_1) \Big \rangle + \Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left|- \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}^Z(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

The wave function $\phi_{1s}^Z(\vec{r}-\vec{r}_1)$ is an eigenfunction of the atomic orbital Hamiltonian in the first term $H\phi_{1s}^Z(\vec{r}-\vec{r}_1) = E_1 \phi_{1s}^Z(\vec{r}-\vec{r}_1)$, so the first term is easily evaluated,

\begin{equation} H_{11}=E_1 - \Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left| \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}^Z(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

The second term has a singularity at $\vec{r}_2$ which makes it difficult to evaluate numerically. We break the second term into an integral over a spherical volume of radius $\delta$ centered around $\vec{r}_2$ and a second integral outside that volume.

\begin{equation} H_{11}=E_1 - \int\limits_{|\vec{r}-\vec{r}_2| < \delta} \phi_{1s}^Z(\vec{r}-\vec{r}_1) \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \phi_{1s}^Z(\vec{r}-\vec{r}_1) d^3r - \int\limits_{|\vec{r}-\vec{r}_2| > \delta} \phi_{1s}^Z(\vec{r}-\vec{r}_1) \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \phi_{1s}^Z(\vec{r}-\vec{r}_1) d^3r. \end{equation}

Close to $\vec{r}_2$, $\exp\left(\frac{ -Z|\vec{r}-\vec{r}_1|}{a_0}\right)\approx \exp\left(\frac{ -Z|\vec{r}_2-\vec{r}_1|}{a_0}\right)$. Using this approximation, the first integral which includes the singularity can be performed analytically for small $\delta$.

\begin{equation} H_{11}=E_1 - \frac{Z^4e^2\delta^2}{2\pi\epsilon_0a_0^3} \exp(-2Z|\vec{r}_1-\vec{r}_2|/a_0) - \int H(|\vec{r}-\vec{r}_2|-\delta ) \phi_{1s}^Z(\vec{r}-\vec{r}_1) \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \phi_{1s}^Z(\vec{r}-\vec{r}_1) d^3r, \end{equation}

The second integral integrates over all space but a Heaviside step function has been introduced. $H(|\vec{r}-\vec{r}_2|-\delta ) = 0$ for $|\vec{r}-\vec{r}_2| < \delta$ and is 1 otherwise. The second integral contains no singularity and can be evaluated numerically.

$H(\|\vec{r}-\vec{r}_2\|-\delta )\frac{a_0^3e\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_1)}{4\pi\epsilon_0 \|\vec{r}-\vec{r}_2\|}$
	$\frac{x}{a_0}$
The integrand of the matrix element plotted along the $x$-axis for $\delta = a_0/10$.

The code below uses a Monte-Carlo method to and calculate $H_{11}$.

<script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = Math.PI; var ec = 1.60217733E-19; //electron charge var eps0 = 8.854187817E-12; //permittivity constant var Z_eff = 1; var delta = a0/(10Z_eff); var E1 = -Z_effZ_eff2.17987e-18; // E1 in Joules function H(x,y,z,x1,y1,z1,delta) { //Heaviside step function r = Math.sqrt(Math.pow(x-x1,2)+Math.pow(y-y1,2)+Math.pow(z-z1,2)); d = r - delta; if (d<0) {hv = 0;} if (d==0) {hv = 0.5;} if (d>0) {hv = 1;} return hv; } 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)/(piMath.pow(a0,3)))Math.exp(-Zar/a0); } document.write("E<sub>1</sub> = "+E1/ec+" eV<br />"); x1 = -0.38E-10; x2 = 0.38E-10; C_int = -(Math.pow(Z_eff,4)ececdeltadelta/(2pia0a0a0eps0))Math.exp(-2Z_eff(x2-x1)/a0); document.write("For Z_effdelta/a0 = "+(Z_effdelta/a0).toPrecision(2)+", the integral over the singularity is, "+C_int/ec+" eV <br />"); L = 8a0/Z_eff; H11 = 0; N = 1000000; //number of random points for (i=0; i<N; i++) { x = L(0.5-Math.random()); y = L(0.5-Math.random()); z = L(0.5-Math.random()); H11 += -H(x,y,z,x2,0,0,delta)ao_1s(x,y,z,x1,0,0,1)ao_1s(x,y,z,x1,0,0,1)/Math.sqrt(Math.pow(x-x2,2)+yy+zz); } H11 = Math.pow(L,3)H11ec/(4pieps0*N); document.write("The numerical integral is = "+H11+" eV<br />"); H11 += C_int/ec + E1/ec; document.write("H<sub>11</sub> = "+H11+" eV"); </script>

Similarly, the off-diagonal Hamiltonian matrix element of a homonuclear diatomic molecule (H₂, O₂, N₂, etc.) with two 1s orbitals located at positions $\vec{r}_1$ and $\vec{r}_2$ is,

\begin{equation} H_{12}=\Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}- \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{1s}^Z(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

This can be broken into two terms,

\begin{equation} H_{12}=\Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left|- \frac{\hbar^2}{2m}\nabla^2 - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|}\right|\phi_{1s}^Z(\vec{r}-\vec{r}_2) \Big \rangle + \Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left|- \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \right|\phi_{1s}^Z(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

The wave function $\phi_{1s}^Z(\vec{r}-\vec{r}_2)$ is an eigenfunction of the atomic orbital Hamiltonian in the first term $H\phi_{1s}^Z(\vec{r}-\vec{r}_2) = E_1 \phi_{1s}^Z(\vec{r}-\vec{r}_2)$, so the first term is easily evaluated,

\begin{equation} H_{12}=E_1S_{12} - \Big \langle \phi_{1s}^Z(\vec{r}-\vec{r}_1)\left| \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \right|\phi_{1s}^Z(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

The second term has a singularity at $\vec{r}_1$ which makes it difficult to evaluate numerically. We break the second term into an integral over a spherical volume of radius $\delta$ centered around $\vec{r}_1$ and a second integral outside that volume.

\begin{equation} H_{12}=E_1S_{12} - \int\limits_{|\vec{r}-\vec{r}_1| < \delta} \phi_{1s}^Z(\vec{r}-\vec{r}_1) \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \phi_{1s}^Z(\vec{r}-\vec{r}_2) d^3r - \int\limits_{|\vec{r}-\vec{r}_1| > \delta} \phi_{1s}^Z(\vec{r}-\vec{r}_1) \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \phi_{1s}^Z(\vec{r}-\vec{r}_2) d^3r. \end{equation}

Close to $\vec{r}_1$, $\exp\left(\frac{ -Z|\vec{r}-\vec{r}_2|}{a_0}\right)\approx \exp\left(\frac{ -Z|\vec{r}_1-\vec{r}_2|}{a_0}\right)$. Using this approximation, the first integral which includes the singularity can be performed analytically for small $\delta$.

\begin{equation} H_{12}=E_1S_{12} - \frac{Z^4e^2\delta^2}{2\pi a_0^3\epsilon_0} \exp(-Z|\vec{r}_1-\vec{r}_2|/a_0) - \int H(|\vec{r}-\vec{r}_1|-\delta ) \phi_{1s}^Z(\vec{r}-\vec{r}_1) \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \phi_{1s}^Z(\vec{r}-\vec{r}_2) d^3r, \end{equation}

The second integral integrates over all space but a Heaviside step function has been introduced. $H(|\vec{r}-\vec{r}_1|-\delta ) = 0$ for $|\vec{r}-\vec{r}_1| < \delta$ and is 1 otherwise. The second integral contains no singularity and can be evaluated numerically.

$H(\|\vec{r}-\vec{r}_1\|-\delta )\frac{a_0^3e\phi_{1s}(\vec{r}-\vec{r}_1)\phi_{1s}(\vec{r}-\vec{r}_2)}{4\pi\epsilon_0 \|\vec{r}-\vec{r}_1\|}$
	$\frac{x}{a_0}$
The integrand of the matrix element plotted along the $x$-axis for $\delta = a_0/10$.

The code below uses a Monte-Carlo method to and calculate $H_{12}$.

<script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = Math.PI; var ec = 1.60217733E-19; //electron charge var eps0 = 8.854187817E-12; //permittivity constant var Z_eff = 1; var delta = a0/(10Z_eff); var S12 = 0.74; //overlap integral var E1 = -Z_effZ_eff2.17987e-18; // E1 in Joules function H(x,y,z,x1,y1,z1,delta) { //Heaviside step function r = Math.sqrt(Math.pow(x-x1,2)+Math.pow(y-y1,2)+Math.pow(z-z1,2)); d = r - delta; if (d<0) {hv = 0;} if (d==0) {hv = 0.5;} if (d>0) {hv = 1;} return hv; } 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)/(piMath.pow(a0,3)))Math.exp(-Zar/a0); } document.write("E<sub>1</sub>S<sub>12</sub> = "+E1S12/ec+" eV<br />"); x1 = -0.38E-10; x2 = 0.38E-10; C_int = -(Math.pow(Z_eff,4)ececdeltadelta/(2pia0a0a0eps0))Math.exp(-Z_eff(x2-x1)/a0); document.write("For Z_effdelta/a0 = "+(Z_effdelta/a0).toPrecision(2)+", the integral over the singularity is, "+C_int/ec+" eV <br />"); H12 = 0; N = 1000000; //number of random points for (i=0; i<N; i++) { x = 8a0(0.5-Math.random())/Z_eff; y = 8a0(0.5-Math.random())/Z_eff; z = 8a0(0.5-Math.random())/Z_eff; H12 += -H(x,y,z,x1,0,0,delta)ao_1s(x,y,z,x1,0,0,Z_eff)ao_1s(x,y,z,x2,0,0,Z_eff)/Math.sqrt(Math.pow(x-x1,2)+yy+zz); } H12 = Math.pow(8a0/Z_eff,3)H12ec/(4pieps0N); document.write("The numerical integral is = "+H12+" eV<br />"); H12 += C_int/ec + S12*E1/ec; document.write("H<sub>12</sub> = "+H12+" eV"); </script>

The elements of the Hamiltonian matrix and overlap matrix were recalcuated in matlab. The scripts they used were ao_1s.m, ao_2s.m, Laplace1s.m, Laplace2s.m, Heavysidestep.m, Hydrogen_1s.m, and Hydrogen_2s.m.
Another Matlab script to calculate the overlap matrix and the Hamiltonian matrix is calcHamiltonian.m. This script uses the atomic orbitals defined in atorbit.m.