This is the project work of Konstantin Rygol for the course PHYS291 - Data Handling in Physics. Hold in Spring semester 2019 at UIB.

The first part of the project will be to calculate the wave functions of the Hydrogen atom. Afterwards several orbitals of the Hydrogen atom will be plotted using the wave functions. For the calculation and graphical representation C++/ROOT are going to be the tools of choice. Due to high performance needs. The plots base on simulations having 10 million grid points.

The Wave function will not be derived here but is assumed to be known by the reader. The electron wave function is given by: \[ \psi_{nlm} = \sqrt{ \left (\frac{2}{na_0} \right)^3 \frac{(n-l-1)!}{2n (n+l)!}} e^{-\rho / 2} \rho^l L^{2l+1}_{n-l-1}(\rho) Y^m_l(\theta,\phi) \]

L is the associated Laguerre polynomial, Y the spherical harmonic and \(a_0\) is the Bohr radius. The spherical coordinates \(\rho\), \(\theta\), \(\phi\) correspond to 2r/n, the polar angle and the azimuthal angle . The 3 quantum numbers are n the principle quantum number, l the angular momentum quantum number and m the magnetic quantum number. The density at a given point is then given by the square of \(\psi\).

The simulation code is written in C++ and can be found in project.cpp . The density has three parameters for implementation reasons a 1D array is chosen with a mapping function that maps the three spherical parameters to the 1Dimensional array. 1 Dimensional arrays offer better control regarding functions and parallelization. The code could be easily parallelized with OPEN-MP. The wave functions are complex but the density is real. C++ is a good choice at it supports complex numbers within its standard library.

The simulation code is built around the following functions:

double getX(double r,double theta ,double phi)

double getY(double r,double theta ,double phi)

double getZ(double r,double theta ,double phi)

indDen(int i,int j , int k,int r_points, int pola_points)

factorial(int n)

std::complex

PSI(double r, double theta, double phi, int n, int l,int m) std::vector

genDensity(double r_max, int r_points, int aze_points, int pola_points,int n,int l, int m) std::vector

calcSurface(std::vector density, double r_max, int r_points, int aze_points, int pola_points)

In the main function the genDensity function is called. It then calls for every grid point the gen PSI to obtain the expected value at that point. Afterwards the calcSurface is called to find a surface of the density. First the equidensity line that encloses 90% of the points is calculated then the points that form the line are plotted. To find the first 90% values all grid points get sorted after their density valuei. Then the value that divides the set to 90:10 is found. This value is called contour_value.

It pays out here that we have chosen polar coordinates. The idea now is to go a long every ray and find the first value that is bigger than the contour_value. Then this point is added to the to be printed points and a marker is set. After Iterating further there will be a point along the ray where the density is lower than the contour_value. This point is then also added to the to be printed list and the marker is removed. So if after a while we might reach a new point with a density greater than the contour_value. A new point point can be added to the list. This allows to find the surface of more complex orbitals.

*The points above the threshold line are added to the surface. For the ISO plot only the first point and the last point of a local maxima above the threshold line are added. The threshold line is at 0.6 in this plot.*

This method is fast as it consists of simple comparisons along a continuous array in the memory. For performance reasons it is therefore important to find a beneficial mapping function from the triple(r,\(\phi\),\(\theta\)) to the continuous 1D memory the data is stored in. This method however comes with a some disadvantages. The density of grid point is not constant meaning that further away from the origin our density of grid points is lower. This results in a coarser resolution in areas further away from the origin.

*This figure shows the density of the grid points for the simulation. It is visible that the density is higher at the center and along the z axis.*

To plot the results two techniques will be used one will be to just plot all points returned by the create surface function the other one will be a a 3d histogram using the option ISO to draw a contour surface. Both techniques have their respective strengths and weaknesses. While the 3d histogram gives a sooth picture of the orbital it also hides information that might be covered under the ISO surface. This is especially true for higher quantum numbers. The scatter plot reveals this information but can be difficult to read. It also appears to change the scatter plot every time it is clicked. I suspect that ROOT might not display all points correctly if there are two many points in the graphic.

A key aspect is also that the following plots are of course 2d projections and highly depend on the perspective. To allow a better comparison all plots are presented here as the standard projection ROOT chooses.

*The 1s orbital look like expected.* *The scatter plot looks it is more squared but this might be a projection mistake.*

*In this plot there seems to be a region around the origin marked where there should be a gap between a lower and an upper half. This can be explained by the projection error.* *The ISO plot shows both a lower and an upper part of the orbital but it is not able to see the difference between the lower and the upper orbital.*

*In this plot there seems to be a region around the origin marked where there should be a gap between a lower and an upper half. This can be explained by the projection error. This plot is similar to the 2p orbitals but in the scatter plot differences can be found.* *The ISO plot shows both a lower and an upper part of the orbital but it is not able to see the difference between the lower and the upper orbital. The ISO plot is not by eye distinguishable from the the 2p orbital.*

*The separation in a lower and upper region is visible . No substructures can be distinguished.*

*The separation in a lower and upper region is visible. It is also visible that there are multiple rings like expected.*

*The plot shows the ringlike structure of the \(4f_3\) orbital. Some substructures can be seen.* *The plot shows the ringlike structure of the \(4f_3\) orbital.*

This was the first bigger project using root as plot program. The satisfaction with root is mixed on the one side there is the outstanding performance over for example matplotlib (a python package) but on the other side the community is small the information found on the internet is sometimes so old that the web page behind a web link has been deleted. The two graphical represenation chosen where the best ones i found in the ROOT documentation. But they are less than ideal.

It is difficult to web search efficiently for ROOT questions. Also other experiences on the web agree that ROOT is more a nique product for people who do have to deal with a lot of data and also might have to plot a lot of data and not every scientist. The integration within C++ code was tried but I failed in getting the right functions from the libraries. Extensive web search did not produce an easy tutorial or solution for that. I therefore used ROOT separately from my simulation code. I for myself conclude that ROOT was a nice experience but will not be part of my future work flow. I am thank full to my teachers and tutors for holding this course.

The code is published under the MIT License. Thanks to killercup for providing the pandoc.css.