The 2p_z atomic orbital is,

\[\begin{equation} \phi^Z_{2p_z}= \sqrt{\frac{Z^3}{\pi 2^5 a_0^3}}\frac{Zz}{a_0} \exp\left(-\frac{Zr}{2a_0}\right). \end{equation}\]

To account for the screening due to the core electrons, the effective nuclear charge $Z$ should be determined by Slater's rules. For carbon, the effective nucelar charge is $Z=3.25$.

$\left(\frac{a_0}{Z}\right)^{3/2}\phi_{2p_z}$
	$\frac{Zz}{a_0}$
The $2p_z$ orbital plotted along the $z$-axis. The maximum is at $z=\frac{2a_0}{Z}$.

To check that the 2p_z 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_{2p_z}=\sqrt{\frac{Z^7}{\pi 2^5 a_0^7}}\left(\frac{Zr}{4a_0}-2\right)\frac{z}{r} \exp\left(-\frac{Zr}{2a_0}\right), \end{equation}\]

and then calculate the energy $E_2= - \frac{\hbar^2}{2m}\frac{\nabla^2\phi_{2p_z}(\vec{r})}{\phi_{2p_z}(\vec{r})} - \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}|}$. The code below defines the 2p_z 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 -35.9 eV at every position for $Z=3.25$.

<script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = 3.1415926535897932; var hbar = 1.05457266E-34; //h-bar var ec = 1.60217733E-19; //electron charge var eps0 = 8.854187817E-12; //permittivity constant var me = 9.1093897E-31; //electron mass var k1 = 0.5*hbar*hbar/me; var k2 = ec*ec/(4*pi*eps0); //Coulomb constant Z_eff = 3.25; 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)); } for (i=0; i<10; i++) { //The energy is calculated for 10 random positions x = 8*a0*(0.5-Math.random())/Z_eff; xi = 8*a0*(0.5-Math.random())/Z_eff; y = 8*a0*(0.5-Math.random())/Z_eff; yi = 8*a0*(0.5-Math.random())/Z_eff; z = 8*a0*(0.5-Math.random())/Z_eff; zi = 8*a0*(0.5-Math.random())/Z_eff; E = -k1*L_ao_2pz(x,y,z,xi,yi,zi,Z_eff)/ao_2pz(x,y,z,xi,yi,zi,Z_eff)-k2*Z_eff/r; document.write("E₂ = "+E/ec+" eV<br>"); } </script>

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

\begin{equation} S_{12}=\langle \phi_{2p_z}^Z(\vec{r}-\vec{r}_1)|\phi_{2p_z}^Z(\vec{r}-\vec{r}_2)\rangle . \end{equation}

For a ethene molecule, $Z_{eff}=3.25$, $\vec{r}_1=-0.67\,\hat{x}$ Å and $\vec{r}_2=0.67\,\hat{x}$ Å. Below $\phi_{2p_z}^Z(\vec{r}-\vec{r}_1)\phi_{2p_z}^Z(\vec{r}-\vec{r}_2)$ is plotted along the $x$-direction at $z=\frac{2a_0}{Z}$ and along the $z$-direction through one of the atoms.

$\left(\frac{a_0}{Z}\right)^3\phi_{2p_z}^Z(\vec{r}-\vec{r}_1)\phi_{2p_z}^Z(\vec{r}-\vec{r}_2)$
	$\frac{Zx}{a_0}\,@\,y=0,\,z=\frac{2a_0}{Z}$

The code below uses a Monte-Carlo method to integrate $\phi_{2p_z}^C(\vec{r}-\vec{r}_1)\phi_{2p_z}^C(\vec{r}-\vec{r}_2)$ and calculate $S_{12}$. For carbon, $Z_{eff}=3.25$.

<script> var a0 = 5.2917721E-11; // Bohr radius in meters var pi = 3.1415926535897932; var Z_eff = 3.25; 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)); } S12 = 0; N = 1000000; //number of random points x1 = -0.67E-10; x2 = 0.67E-10; for (i=0; i<N; i++) { x = 30*a0*(0.5-Math.random())/Z_eff; // -15 to +15 Zr/a0 y = 30*a0*(0.5-Math.random())/Z_eff; z = 30*a0*(0.5-Math.random())/Z_eff; S12 += ao_2pz(x,y,z,x1,0,0,Z_eff)*ao_2pz(x,y,z,x2,0,0,Z_eff); } document.write("S<sub>12<\/sub> = "+Math.pow(30*a0/Z_eff,3)*S12/N); </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_{2p_z}(\vec{r})|\phi_{2p_z}(\vec{r})\rangle$ which should equal to 1 if the wave function is properly normalized.

The Hamiltonian matrix element $H_{11}$ for ethene with two 2p_z orbitals is,

\begin{equation} H_{11}=\Big \langle \phi_{2p_z}^C(\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_{2p_z}^C(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

This can be broken into two terms,

\begin{equation} H_{11}=\Big \langle \phi_{2p_z}^C(\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_{2p_z}^C(\vec{r}-\vec{r}_1) \Big \rangle + \Big \langle \phi_{2p_z}^C(\vec{r}-\vec{r}_1)\left|- \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_2|} \right|\phi_{2p_z}^C(\vec{r}-\vec{r}_1) \Big \rangle . \end{equation}

The wave function $\phi_{2p_z}^C(\vec{r}-\vec{r}_1)$ is an eigenfunction of the atomic orbital Hamiltonian in the first term $H\phi_{2p_z}^C(\vec{r}-\vec{r}_1) = E_2 \phi_{2p_z}^C(\vec{r}-\vec{r}_2)$, so the first term is easily evaluated,

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

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 = 3.1415926535897932; var ec = 1.60217733E-19; //electron charge var eps0 = 8.854187817E-12; //permittivity constant var E2 = -35.92*ec; // E2 in Joules var Z_eff = 3.25; 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)); } document.write("E<sub>2</sub> = "+E2/ec+" eV<br />"); x1 = -0.67E-10; x2 = 0.67E-10; H11 = 0; N = 1000000; //number of random points for (i=0; i<N; i++) { x = 30*a0*(0.5-Math.random())/Z_eff; // -15 to +15 Zr/a0 y = 30*a0*(0.5-Math.random())/Z_eff; z = 30*a0*(0.5-Math.random())/Z_eff; H11 += -ao_2pz(x,y,z,x1,0,0,Z_eff)*ao_2pz(x,y,z,x1,0,0,Z_eff)/Math.sqrt(Math.pow(x-x2,2)+y*y+z*z); } H11 = Math.pow(30*a0/Z_eff,3)*H11*ec*Z_eff/(4*pi*eps0*N); document.write("The numerical integral is = "+H11+" eV<br />"); H11 += E2/ec; document.write("H<sub>11</sub> = "+H11+" eV"); </script>

The Hamiltonian matrix element $H_{12}$ for ethene with two 2p_z orbitals located at positions $\vec{r}_1$ and $\vec{r}_2$ is,

\begin{equation} H_{12}=\Big \langle \phi_{2p_z}^C(\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_{2p_z}^C(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

This can be broken into two terms,

\begin{equation} H_{12}=\Big \langle \phi_{2p_z}^C(\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_{2p_z}^C(\vec{r}-\vec{r}_2) \Big \rangle + \Big \langle \phi_{2p_z}^C(\vec{r}-\vec{r}_1)\left|- \frac{Ze^2}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|} \right|\phi_{2p_z}^C(\vec{r}-\vec{r}_2) \Big \rangle . \end{equation}

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

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

Below is a plot of the integrand of the second term. Even though there is a factor of $\frac{1}{|\vec{r}-\vec{r}_1|}$, there is no singularity because the $\phi_{2p_z}^C(\vec{r}-\vec{r}_1)$ wave function goes to zero at that point.

$-\frac{a_0^3}{Z^3}\frac{\phi_{2p_z}^Z(\vec{r}-\vec{r}_1)\phi_{2p_z}^Z(\vec{r}-\vec{r}_2)}{4\pi\epsilon_0 |\vec{r}-\vec{r}_1|}$
	$\frac{Zz}{a_0}$
The integrand of the matrix element plotted along the $z$-direction through the atom at $\vec{r}_1$.

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 = 3.1415926535897932; var ec = 1.60217733E-19; //electron charge var eps0 = 8.854187817E-12; //permittivity constant var S12 = 0.27; //overlap integral var E2 = -35.92*ec; // E2 in Joules var Z_eff = 3.25; 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)); } document.write("E<sub>2</sub>S<sub>12</sub> = "+E2*S12/ec+" eV<br />"); x1 = -0.67E-10; x2 = 0.67E-10; H12 = 0; N = 1000000; //number of random points for (i=0; i<N; i++) { x = 30*a0*(0.5-Math.random())/Z_eff; // -15 to +15 Zr/a0 y = 30*a0*(0.5-Math.random())/Z_eff; z = 30*a0*(0.5-Math.random())/Z_eff; H12 += -ao_2pz(x,y,z,x1,0,0,Z_eff)*ao_2pz(x,y,z,x2,0,0,Z_eff)/Math.sqrt(Math.pow(x-x1,2)+y*y+z*z); } H12 = Math.pow(30*a0/Z_eff,3)*H12*ec*Z_eff/(4*pi*eps0*N); document.write("The numerical integral is = "+H12+" eV<br />"); H12 += S12*E2/ec; document.write("H<sub>12</sub> = "+H12+" eV"); </script>

A Matlab script to calculate the overlap matrix and the Hamiltonian matrix is calcHamiltonian.m. This script uses the atomic orbitals defined in atorbit.m.