Advanced Chemistry - Mono Mole

December 7, 2022December 2, 2024

Analytical solution of the energy equation for He using the Hartree self-consistent field method

The analytical solution of the energy equation of He involves finding analytical expressions for all terms in eq8 using trial one-electron wavefunctions. As per the numerical method, guess values are used for the variables in all trial one-electron wavefunctions except one, e.g. $\phi_1(\boldsymbol{\mathit{r}}_1)$ . Eq8 is then minimised to obtain the solutions for $E$ and $\phi_1(\boldsymbol{\mathit{r}}_1)$ . The process is repeated until $E$ and all $\phi_1(\boldsymbol{\mathit{r}}_1)$ become invariant. For He, eq8 is

$E=H_1+H_2+J_{12}\;\;\;\;\;\;\;\;(48)$

where $H_i=\int\phi_i^*(\mathbf r_i)\left(-\frac{1}{2}\nabla_i^2-\frac{Z}{r_i}\right)\phi_i(\mathbf r_i)d\mathbf r_i$ and $J_{ij}=\int\int\frac{\vert\phi_i\vert^2\vert\phi_j\vert^2}{r_{ij}}d\mathbf r_id\mathbf r_j$ .

Our computation employs the following assumptions:

We base our iteration on the unrestricted case, where electrons $j=1$ and $i=2$ are distinctively expressed by $\phi_1$ and $\phi_2$ respectively.
All terms on the RHS of eq48 are determined using Slater-type orbitals, where

$\phi_1=\sqrt\frac{\alpha^3}{\pi}e^{-\alpha r_1}\;\;\;and\;\;\;\phi_2=\sqrt\frac{\beta^3}{\pi}e^{-\beta r_2}$

Substituting $\phi_1$ in $H_1$ , using the identity $\int_0^{\infty}x^ne^{-mx}dx=\frac{n!}{m^{n+1}}$ for $\left\langle\phi_1\bigg| \small{-\frac{1}{2}\nabla_1^2}\bigg|\phi_1\right\rangle$ and integrating by parts for $\left\langle\phi_1\bigg| \small{\frac{2}{r_1}}\bigg|\phi_1\right\rangle$ , we have

$H_1=\frac{\alpha^2}{2}-2\alpha \;\;\;\;\;\;\;\;(50)$

Similarly, substituting $\phi_2$ in $H_2$ , we have

$H_2=\frac{\beta^2}{2}-2\beta \;\;\;\;\;\;\;\;(51)$

Substituting $\phi_1$ , $\phi_2$ and eq47 in $J_{12}$ and multiplying the resultant equation by $4\pi[Y_0^0(\theta_1,\phi_1)]^*Y_0^0(\theta_2,\phi_2)=1$ , where $Y_0^0(\theta,\phi)=\frac{1}{\sqrt{4\pi}}$ is the spherical harmonics for $l=0$ , we have

$J_{12}=\frac{4\alpha^3\beta^3}{\pi}\int\int e^{-2(\alpha r_1+\beta r_2)}\sum_{n=0}^{\infty}\sum_{m=-n}^n\frac{4\pi}{2n+1}\frac{r_<^n}{r_>^{n+1}}$

$\times Y_n^{\vert m\vert}(\theta_1,\phi_1)[Y_n^{\vert m\vert}(\theta_2,\phi_2)]^*[Y_0^0(\theta_1,\phi_1)]^*Y_0^0(\theta_2,\phi_2)d\boldsymbol{\mathit{r}}_1d\boldsymbol{\mathit{r}}_2$

Due to the orthogonality of the spherical harmonics, the only integral that survives upon expanding the summation is when $n=0$ and $m=0$ . Since, $\int d\boldsymbol{\mathit{r}}=\int\int\int r^2drsin\theta d\theta d\phi$ , the above equation becomes:

$J_{12}=16\alpha^3\beta^3\int_0^{\infty} e^{-2\alpha r_1}r_1^2dr_1\frac{1}{r_>}\int_0^{\infty}e^{-2\beta r_2}r_2^2dr_2$

We proceed by integrating with respect to $r_2$ first. Since $r_2$ ranges from 0 to $\infty$ , we can split the integral into two parts, one from 0 to $x$ and the other from $x$ to $\infty$ . Supposing $x=r_1$ , the above equation becomes:

$J_{12}=16\alpha^3\beta^3\int_0^{\infty} r_1^2e^{-2\alpha r_1}\frac{1}{r_>}\left(\int_0^{r_1}r_2^2e^{-2\beta r_2}dr_2+\int_{r_1}^{\infty}r_2^2e^{-2\beta r_2}r_2^2dr_2\right)dr_1$

As $r_2$ increases from 0 and approaches $r_1$ , it must be less than $r_1$ and so $r_>=r_1$ . Similarly, as $r_2$ increases from $r_1$ to $\infty$ , it must be greater than $r_1$ and so $r_>=r_2$ . Therefore,

$J_{12}=16\alpha^3\beta^3\int_0^{\infty} r_1^2e^{-2\alpha r_1}\bigg(\frac{1}{r_1}\int_0^{r_1}r_2^2e^{-2\beta r_2}dr_2$

$+\int_{r_1}^{\infty}r_2e^{-2\beta r_2}r_2^2dr_2\bigg)dr_1\;\;\;\;\;\;\;\;(52)$

Using the identity $\int_0^b x^n e^{-ax}dx=\frac{n!}{a^{n+1}}\left[1-e^{-ab}\sum_{i=0}^n\frac{(ab)^i}{i!}\right]$ for the first integral within parentheses, integrating by parts for the second integral within parentheses, and employing the identity $\int_0^{\infty} x^n e^{-mx}dx=\frac{n!}{m^{n+1}}$ for the integral with respect to $r_1$ yields

$J_{12}=\frac{\alpha\beta(\alpha^2+3\alpha\beta+\beta^2)}{(\alpha+\beta)^3}\;\;\;\;\;\;\;\;55$

Substituting eq50, eq51 and eq55 in eq48 gives

$E=\frac{\alpha^2}{2}-2\alpha+\frac{\beta^2}{2}-2\beta+\frac{\alpha\beta(\alpha^2+3\alpha\beta+\beta^2)}{(\alpha+\beta)^3}\;\;\;\;\;\;\;\;56$

Differentiating eq56 with respect to $\alpha$ and $\beta$ results in

$\frac{\partial E}{\partial\alpha}=\alpha-2+\frac{4\alpha\beta^3+\beta^4}{(\alpha+\beta)^4}\;\;\;\;\;\;\;\;57$

$\frac{\partial E}{\partial\beta}=\beta-2+\frac{\alpha^4+4\alpha^3\beta}{(\alpha+\beta)^4}\;\;\;\;\;\;\;\;58$

Eq57 and eq58 are used in an iterative algorithm to find $E$ , $\phi_1$ and $\phi_2$ . The procedure is as follows:

Substitute an initial guess value of $\beta$ , e.g., $\beta=2$ , in eq57 and solve for $\alpha$ by setting eq57 equal to zero. We then substitute the solution of $\alpha$ and the initial guess value of $\beta$ in eq56 to find $E$ .
To obtain an improved estimate of $\beta$ , we substitute the solution of $\alpha$ from the previous step in eq58, and solve for $\beta$ by setting eq58 equal to zero. We then substitute the solution of $\beta$ and the value of $\alpha$ found in the previous step in eq56 to find a better estimate of $E$ .
Steps 1 and 2, which form an iteration set, are repeated until the values of $\alpha$ and $\beta$ are invariant up to six decimal points.

Alternatively, we can set up an iterative table in Excel as follows:

Cell C2 is the initial guess value of $\beta$ . The formula for D2:D7 is eq56, i.e. =(B2^2)/2-2*B2+(C2^2)/2-2*C2+(B2*C2*((B2^2)+3*B2*C2+(C2^2)))/((B2+C2)^3). The formulae for B3, C4, B5, C6 and B7 are =B2, =C3, =B4, =C5 and =B6, respectively. To compute the first iteration set, we employ the Excel Solver application to minimise D2 with respect to B2, with the following settings:

Set objective: D2
To: Min
By changing variable cells: B2
Selecting a solving method: GRG Nonlinear

We then solve for D3 by changing the ‘Set objective’ field and ‘By changing variable cells’ field to D3 and C3 respectively. This procedure is repeated for subsequent iteration sets until the values of $\alpha$ and $\beta$ are invariant up to six decimal points.

The results are as follows:

The final theoretical value of $E=-2.8476562$ has a deviation of about 1.93% versus the experiment data of the ground state of helium.

Next article: Antisymmetriser

Previous article: The expansion of $\frac{1}{r_{ij}}$ (the reciprocal of the inter-electron distance)

Content page of computational chemistry

Content page of advanced chemistry

Main content page

December 7, 2022September 26, 2024

Analytical solution of the energy equation for Li using the Hartree-Fock method

The analytical solution of the energy equation for Li involves finding analytical expressions for all terms in eq94 using one-electron trial wavefunctions. Guess values are used for the variables in all the wavefunctions except one, e.g. $\phi_1(\boldsymbol{\mathit{ r}}_1)$ . Eq94 is then minimised to obtain the solutions for E and $\phi_1(\boldsymbol{\mathit{ r}}_1)$ . The process is repeated until E and all $\phi_i(\boldsymbol{\mathit{ r}}_i)$ become invariant.

After removing $K_{12}$ and $K_{23}$ due to spin orthogonality from eq94, we have

$E=H_1+H_2+H_3+J_{12}+J_{13}+J_{23}-K_{13}\;\;\;\;\;\;\;\;(113)$

where
$H_i=\int\phi_i^*(\boldsymbol{\mathit{r}}_i)\left(-\frac{1}{2}\nabla_i^2-\frac{Z}{r_i}\right)\phi_i(\boldsymbol{\mathit{r}}_i)d\boldsymbol{\mathit{r}}_i$
$J_{ij}=\int\int\phi_i^*(\boldsymbol{\mathit{r}}_i)\phi_j^*(\boldsymbol{\mathit{r}}_j)\frac{\phi_i(\boldsymbol{\mathit{r}}_i)\phi_j(\boldsymbol{\mathit{r}}_j)}{r_{ij}}d\boldsymbol{\mathit{r}}_id\boldsymbol{\mathit{r}}_j$
$K_{ij}=\int\int\phi_i^*(\boldsymbol{\mathit{r}}_i)\phi_j^*(\boldsymbol{\mathit{r}}_j)\frac{\phi_j(\boldsymbol{\mathit{r}}_i)\phi_i(\boldsymbol{\mathit{r}}_j)}{r_{ij}}d\boldsymbol{\mathit{r}}_id\boldsymbol{\mathit{r}}_j$

Our computation employs the following assumptions:

The ground state of Li is described by a single Slater determinant.
Eq113 is minimised iteratively via the restricted open-shell method, where the 1s orbital is doubly occupied with $\phi_1=\phi_2=\sqrt{\frac{\zeta_1^3}{\pi}}e^{-\zeta_1r}$ , while the 2s orbital is singly occupied. Therefore, we can substitute eq56 in eq113, where $\alpha=\beta=\zeta_1$ , to give

$E=\zeta_1^2-\frac{27}{8}\zeta_1+H_3+J_{13}+J_{23}-K_{13}\;\;\;\;\;\;\;\;(114)$

The 2s Slater-type orbital $\phi_3^{' }$ is derived via the Gram-Schmidt process, as the generic 2s Slater-type orbital is not orthogonal to the 1s Slater-type orbital.

Question

Using the normalised 1s Slater-type orbital $\phi_1=\sqrt{\frac{\zeta_1^3}{\pi}}e^{-\zeta_1r}$ and the un-normalised generic 2s Slater-type orbital $\chi_3=re^{-\zeta_3r}$ , derive $\phi_3^{' }$ via the Gram-Schmidt process.

Answer

With reference to eq111, we have

$\phi_3^{' }=N'\left(\chi_3-\frac{\int\chi_3\phi_1d\tau}{\int\phi_1\phi_1d\tau}\phi_1\right)=N'\left(\chi_3-\int\chi_3\phi_1d\tau\phi_1\right)$

where N’ is the normalisation constant.

Substituting $\phi_1$ and $\chi_3$ in the above equation and integrating using the identity $\int_0^{\infty}x^ne^{-mx}dx=\frac{n!}{m^{n+1}}$ , and then normalising and rearranging the result, we have

$\phi_3^{' }=N\phi_3-M\phi_1\;\;\;\;\;\;\;\;(116)$

where
$N=\sqrt{\frac{(\zeta_1+\zeta_3)^8}{(\zeta_1+\zeta_3)^8-192\zeta_1^3\zeta_3^5}}$
$\phi_3=\sqrt{\frac{\zeta_3^5}{3\pi}re^{-\zeta_3r}$
$M=\frac{8}{(\zeta_1+\zeta_3)^4}\sqrt{\frac{3\zeta_1^3\zeta_3^5(\zeta_1+\zeta_3)^8}{(\zeta_1+\zeta_3)^8-192\zeta_1^3\zeta_3^5}}$

Let’s evaluate the remaining terms in eq114, beginning with $H_3=\int\phi_3^{'*}\left(-\frac{1}{2}\nabla^2-\frac{Z}{r}\right)\phi_3^{'}d\tau$ . Substituting eq116 in $H_3$ and integrating,

$H_3=N^2\left(\frac{\zeta_3^2}{6}-\frac{\zeta_3Z}{2}\right)+NM\frac{\sqrt{48\zeta_1^7\zeta_3^5}}{(\zeta_1+\zeta_3)^4}-\frac{NM}{(\zeta_1+\zeta_3)^3}\sqrt{\frac{64\zeta_1^5\zeta_3^5}{3}}$ $+\frac{NMZ}{(\zeta_1+\zeta_3)^3}\sqrt{\frac{256\zeta_1^3\zeta_3^5}{3}}+NM\frac{\sqrt{48\zeta_1^3\zeta_3^9}}{(\zeta_1+\zeta_3)^4}-\frac{NM}{(\zeta_1+\zeta_3)^3}\sqrt{\frac{256\zeta_1^3\zeta_3^7}{3}}$ $+\frac{NM}{(\zeta_1+\zeta_3)^2}\sqrt{\frac{16\zeta_1^3\zeta_3^5}{3}}+M^2\frac{\zeta_1^2}{2}-M^2Z\zeta_1\;\;\;\;\;\;\;\;(119)$

Next, we evaluate $J_{13}$ , which is equal to $J_{23}$ . Substituting eq116 and eq47 in $J_{13}=\int_0^{\infty}\phi_1^*(\boldsymbol{\mathit{r}}_1)\phi_3^{'*}(\boldsymbol{\mathit{r}}_3)\frac{\phi_1(\boldsymbol{\mathit{r}}_1)\phi_3^{' }(\boldsymbol{\mathit{r}}_3)}{r_{13}}d\tau$ and integrating (refer to the integration procedure in this article), we have

$J_{13}=J_{23}=\frac{N^2\zeta_3}{2}-\frac{N^2\zeta_3^5}{2(\zeta_1+\zeta_3)^4}-\frac{N^2\zeta_1\zeta_3^5}{(\zeta_1+\zeta_3)^5}-\frac{2NM}{(\zeta_1+\zeta_3)^3}\sqrt{\frac{64\zeta_1^3\zeta_3^5}{3}}$ $+\frac{2NM}{(3\zeta_1+\zeta_3)^3}\sqrt{\frac{64\zeta_1^3\zeta_3^5}{3}}+\frac{NM\sqrt{768\zeta_1^5\zeta_3^5}}{(3\zeta_1+\zeta_3)^4}+\frac{5}{8}M^2\zeta_1\;\;\;\;\;\;\;\;(129)$

Finally, we evaluate $K_{13}=\int_0^{\infty}\phi_1^*(\boldsymbol{\mathit{r}}_1)\phi_3^{'*}(\boldsymbol{\mathit{r}}_3)\frac{\phi_1(\boldsymbol{\mathit{r}}_3)\phi_3^{ '}(\boldsymbol{\mathit{r}}_1)}{r_{13}}d\tau$ . Referencing the steps in integrating $J_{13}$ , we have

$K_{13}=\frac{44N^2\zeta_1^3\zeta_3^5}{(\zeta_1+\zeta_3)^7}-\frac{NM}{(\zeta_1+\zeta_3)^3}\sqrt{\frac{64\zeta_1^3\zeta_3^5}{3}}+\frac{NM}{(3\zeta_1+\zeta_3)^3}\sqrt{\frac{64\zeta_1^3\zeta_3^5}{3}}$ $+NM\frac{\sqrt{192\zeta_1^5\zeta_3^5}}{(3\zeta_1+\zeta_3)^4}-NM\frac{\sqrt{192\zeta_1^5\zeta_3^5}}{(\zeta_1+\zeta_3)^4}+\frac{NM}{(\zeta_1+\zeta_3)^4}\frac{\sqrt{3072\zeta_1^9\zeta_3^5}}{(3\zeta_1+\zeta_3)^2}$ $+\frac{NM}{\sqrt 3(\zeta_1+\zeta_3)^3}\frac{\sqrt{16384\zeta_1^9\zeta_3^5}}{(3\zeta_1+\zeta_3)^3}+\frac{NM}{(\zeta_1+\zeta_3)^2}\frac{\sqrt{3072\zeta_1^9\zeta_3^5}}{(3\zeta_1+\zeta_3)^4}+\frac{5}{8}M^2\zeta_1\;\;\;\;\;\;\;\;(142)$

Substituting eq119, eq129, eq142 and the atomic number $Z=3$ in eq114, we have the analytical expression of E, which is used in an iterative algorithm to determine the ground state of Li. The procedure involves setting up an iterative table in Excel as follows:

Cell C2 is the initial guess value of $\zeta_3$ . Cell F2 contains the entire formula of the analytical expression of E, with links to B1, C1, D1 and E1. Cell D2 and E2 contain the formulae of N and M respectively, with links to B1 and C1. The formulae for B3, D3 and E3 are ‘=B2’, ‘=D2’ and ‘=E2’ respectively. To compute the first iteration set, we employ the Excel Solver application to minimise F2 with respect to B2, with the following settings:

Set objective: F2
To: Min
By changing variable cells: B2
Selecting a solving method: GRG Nonlinear

We then solve for F3 by changing the ‘Set objective’ field and ‘By changing variable cells’ field to D3 and C3 respectively. This procedure is repeated for subsequent iteration sets until the values of $\zeta_1$ and $\zeta_3$ are invariant up to seven decimal points.

The results are as follows:

The final theoretical value of $E=-7.4179188$ has a deviation of about 0.80% versus the experiment data of the ground state of Li.

Next article: Analytical solution of the Hartree-Fock method for Be, B and C

Previous article: Exchange integral