Calculate an ionized hydrogen molecule¶

In this tutorial, we will calculate the properties of the dihydrogen cation ($ H_2^+ $). It is the simplest molecular ion, as it consists of two hydrogen nuclei and one electron. The simplicity makes it an ideal model system and is the reason why it has been extensively studied.

Exploring 3D Calculations¶

To analyze the properties of $ H_2^+ $, we must solve a 3D quantum mechanical problem. Let’s begin with a simple example calculation in 3D to understand how the computational library handles this scenario and to observe the differences in approach compared to simpler systems.

When transitioning from 1D or 2D to 3D calculations, the primary change is the parameter sdims=3, which specifies the calculation in three spatial dimensions. In fact, this is the only modification required. You can experiment with this in our live environment. Simply input the code, adjust sdims to 1 or 2, and the program will execute the calculations for the chosen dimensions.

In 3D, the output of sg.quickplot changes significantly. It now displays isosurfaces of the wave function of different probabilities. You can toggle the visibility of specific probability levels by a click on the corresponding label. One thing to note: for the given potential, the grid size may be slightly too small, causing the wave function to not fully decay to zero at the boundaries. This results in inaccurate energy values, but for now, the focus is on demonstrating the plotting capability.

Calculating the Dihydrogen Cation Properties¶

To calculate the properties of the dihydrogen cation ($H^+_2$), we first need to define a Hamiltonian that accurately describes the system. This is done using the Born–Oppenheimer approximation, which provides a non-relativistic, time-independent molecular Hamiltonian. The Hamiltonian includes the energy contributions from electrons ($e^-$), nuclei ($n^+$), and their interactions:

\begin{equation} {\cal H} = - \underbrace{ \sum_i {\frac{1}{2} \nabla^2_i} }_{e^- \textrm{ kinetic energy }} - \underbrace{ \sum_A {\frac{1}{2M_A} \nabla^2_A} }_{n^+ \textrm{ kinetic energy}} - \underbrace{ \sum{\frac{Z_A}{r_{i,A}}} }_{e^- n^+\textrm{ attraction}} + \underbrace{ \sum_{i>j}{\frac{1}{r_{ij}}} }_{e^-\textrm{ repulsion}} + \underbrace{ \sum_{B>A} \frac{Z_A \cdot Z_B}{Z_{AB}} }_{n^+\textrm{ repulsion}} \end{equation}

• Electron Kinetic Energy: This is the energy associated with the motion of the electron, where $\nabla^2_i$ represents the second derivative in each spatial dimension.
• Nucleus Kinetic Energy: Similar to the electron, the second term describes the kinetic energy of the nucleus. It includes the mass of the nucleus which is significantly larger than the electron mass ($M_A$ >> $M_e$). Thus, it is often neglected, except when calculating nuclear vibrational modes.
• Electron/Nuclues Attraction: This term accounts for the attractive Coulomb force between the electron and the nucleus. Here, $Z_A$ is the nuclear charge, and $r_{i,A}$ is the distance between the electron $i$ and the nucleus $A$.
• electron/electron repulsion: This energy results from the repulsive Coulomb interaction between electrons, where $r_{i,j}$ is the distance between the two electrons.
• nucleus/nucleus repulsion: This term represents the repulsive force between the nuclei due to their positive charges.

Simplifying the Hamiltonian for $H^+_2$¶

For $H^+_2$, the Hamiltonian simplifies considerably:

• With only one electron, the first term reduces to the Laplacian operator (sum of second derivatives in all three dimensions).
• The nucleus kinetic energy is neglected, setting that term to zero.
• The electron-nucleus attraction includes contributions from both nuclei (A and B).
• There is no electron-electron repulsion because only one electron is present.
• Finally, the nucleus-nucleus repulsion remains as a significant term.

Thus, the Hamiltonian for $H^+_2$ becomes:

\begin{equation} {\cal H_{h2+}} = -\frac{1}{2} \Delta_i - \frac{Z_A}{r_{i,A}} - \frac{Z_B}{r_{i,B}} + \frac{Z_A \cdot Z_B}{Z_{AB}} \end{equation}

Fortunately, these terms do not need to be manually implemented. SigSpace automatically accounts for these contributions when atoms are added to the system. Understanding these terms, however, is crucial for interpreting the results.

Solving a 3D Problem for Two Hydrogen Atoms¶

Now, let’s solve a 3D problem involving two hydrogen atoms to compute the properties of $H^+_2$. This involves using the ${\cal H_{h2+}}$ Hamiltonian and leveraging the library’s built-in capabilities to efficiently calculate and visualize the system's wave function and energy levels.

Let's walk through the code:

• We define the distances between the hydrogen atoms and store it in L_list. In this example we use 20 points logarithmically spaced between 0.7 and 8. The choice of logspace (instead of linspace) allows for a denser sampling at shorter distances, where the results change more rapidly, as we will observe shortly.

• We create a 3D Hamiltonian using H = sg.Hamiltonian(sdims=3) and add two hydrogen atoms with sg.Atom. The atoms are positioned along the x-axis, separated by a distance L. As described earlier, the solver automatically calculates the potential based on the Born-Oppenheimer approximation for each atom added.

• After solving the problem, the calculated energy is stored in an array: E_list.append(float(sol.energy)). Note that sol.energy is a Metric object, which includes the unit of the value. To access only the numerical value, it is cast to a float.

Finally, we plot the calculated energies as a function of the distance between the hydrogen atoms:

Understanding the Energy Curve¶

The characteristic shape in the figure above arises from the contributions in the Hamiltonian. The potential in ${\cal H_{h2+}}$ includes two attractive terms (electron-nucleus attractions) and one repulsive term (nucleus-nucleus repulsion).

At very small distances between the nuclei, the energy increases rapidly due to the strong repulsive force between the positively charged nuclei.

As the distance increases, the energy reaches a minimum around $L=2$. In this region the electron-nucleus attraction becomes most significant, resulting in a minimum energy. Since the system naturally favors the lowest energy configuration, this minimum is the bond length of the molecule.

At larger distances the energy increases again and converges to the value for a single hydrogen atom, as the two nuclei are effectively isolated. The energy therefore saturates at $-0.5$ Ha, which is the energy of a single hydrogen atom (see the tutorial on atomic systems for more details).

Thus, we can extract the optimal bond length from the curve:

The energy at the same location on the curve is called the ionization energy. It describes the amount of energy needed to get a free electron.

You can also calculate the dissociation energy of the dihydrogen cation ($H^+_2$). This energy represents the amount required to break the bond between the two hydrogen nuclei. It is determined by the difference between the minimum energy (at the optimal bond length) and the energy of the system when the nuclei are infinitely separated, the free particle case: