## Harmonic Oscillator Basis Set
<br>
Solving the Schrödinger equation with a potential energy function defined as a quantum Harmonic Oscillator yields a set of wavefunctions that can be used as basis set for quantum mechanical calculations. The Schrödinger equation with the Harmonic Oscillator is defined as follows: 
<br><br>
$$\begin{align}
    \hat{H}\Psi &= E\Psi &&\mbox{Schrödinger Equation}\\
    \\
    \hat{H} &= \hat{T} + \hat{V} &&\mbox{Hamiltonian Operator Expansion}\\
    \\
    \hat{T} &= -\frac{\hbar^2}{2m}\nabla^2 &&\mbox{Kinetic Energy Operator} \\
    \hat{T} &= -\frac{\hbar^2}{2m} \frac{\partial^2}{dr^2} &&\mbox{Simplifying KE Operator to 1D space} \\
    \\
    \hat{V} &= \frac{1}{2}kr^2 &&\mbox{Potential Energy Operator = Harmonic Oscillator Potential} \\
    \hat{H} &= -\frac{\hbar^2}{2m}\frac{\partial^2}{\partial r^2} + \frac{1}{2}kr^2
\end{align}$$
<br>
The wavefunctions produced are:
<br><br>
$$\begin{align}
    \psi_0 &= c_0 \cdot e^{-\frac{\alpha x^2}{2}} \\
    \psi_1 &= c_0 \cdot \sqrt{2\alpha}y \cdot e^{-\frac{\alpha x^2}{2}} \\ 
    \psi_2 &= c_0 \cdot \frac{2y^2 - 1}{\sqrt{2}} \cdot e^{-\frac{\alpha x^2}{2}} \\
    \psi_3 &= c_0 \cdot \frac{2y^3 - 3y}{\sqrt{3}} \cdot e^{-\frac{\alpha x^2}{2}} \\
    \vdots \\
    \psi_n &= c_0 \cdot H_n(y) \cdot e^{-\frac{\alpha x^2}{2}} \\
    \\
    c_0 &= \left[\frac{\alpha}{\pi}\right]^\frac{1}{4} \\
    y &= \sqrt{\alpha}x \\
    x &= r - r_e \\
    r &= \mbox{Bond Distance} \\
    r_e &= \mbox{Optimal Internuclear Bond Distance} \\
    \alpha &= \frac{\mu\omega}{\hbar} \\
    \mu &= \hbox{Reduced Atomic Mass} \\
    \omega &= \hbox{First Term Vibrational Constant} \\
    H_n(y) &= \hbox{The nth Hermite Polynomial} \\
    \end{align}$$
<br>
The Hermite polynomials are added into the wavefunction to allow $\psi$ to be orthonormal and satisfy the orthonormality condition stated below: 
$$\int_{-\infty}^{\infty}{\psi_i \psi_j} = \delta_{i,j}$$
<br>
A sample of Hermite polynomials for selected $n$ values are presented below:
$$\begin{matrix}
    n & H_n \\
    0 & 1 \\
    1 & 2y \\
    2 & 4y^2-2 \\
    3 & 8y^3 - 12y 
\end{matrix}$$
<br>
All information in this notebook originates from the following <a href="http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hosc5.html">website</a>, and the <i>Franck-Condon Calculations</i> Mathematica notebook provided by Dr. Jerry LaRue of Chapman University.    


```python
#Write out terms for KE, and PE
#V = 0 Wavefunction
#Move to other notebook and make like a worksheet
```

## Unit Analysis
<br>

The question of units used for the wavefunction along with its inputs are of the utmost importantance. Many of the more confusing details of the wavefunction implementation are related towards maintaining the correct units for the wavefunction. <br>
To begin with, the overlap integral of two wavefunctions: $\int_{-\infty}^{\infty}{\psi_i\psi_j}=\delta_{ij}$, should result in a unitless quantity of 1 or 0. Unitless normalization of the function can handle ensuring that the output of this integral can only be 0 or 1, but the units for the parameters of the wavefunction must be chosen carefully to allow the integral to be unitless. <br>
The integral will result in units of $r\cdot\psi^2$, normally $r$ will have units of $\mathring A$. Thus, it would imply that $\psi^2$ must have units of $\frac{1}{\mathring A}$, further implying that $\psi$ will have units of $\frac{1}{\sqrt{\mathring A}}$ as: 
$$
\begin{align}
    r\cdot &\psi^2 \\
    \mathring A\cdot &\frac{1}{\sqrt{\mathring A}}^2 \\
    \mathring A\cdot &\frac{1}{\mathring A} \\
    &\frac{\mathring A}{\mathring A} \\
    \mbox{Unitless }&\mbox{Quantity}
\end{align}
$$
<br>
One method to ensure that the wavefunction will have units of $\frac{1}{\sqrt{\mathring A}}$ is to have the wavefunction be unitless by creating a term that can be multiplied by the quantities having units in the wavefunction equation to cancel out the units. Then, a unit quantity of $\frac{1}{\sqrt{\mathring A}}$ is multiplied by the unitless output of the wavefunction.<br>
Given this unit scheme, the inputs and variables of the wavefunction should have the following units: 
$$
\begin{align}
    r &= \mathring A \\
    r_e &= \mathring A \\
    x &= \mathring A \\
    \mu &= \mbox{AMU} \\
    \\
    \omega &= \frac{1}{s} &&\mbox{Needs to be converted from }\frac{1}{cm} \\ 
    \\
    \omega &= \frac{1}{cm} \\
    \omega &= \frac{1}{cm} \cdot 100 \\
    \omega &= \frac{1}{m} \cdot c\left[\mbox{Speed of Light}\right] \\
    \omega &= \frac{1}{m} \cdot \frac{m}{s} \\
    \omega &= \frac{1}{s} \\
\end{align}
$$
<br>
$\hbar$ normally has units of $\frac{kg m^2}{s}$, but in order to properly calculate $\alpha$, the units of $\hbar$ will be converted to $\frac{AMU \mathring A^2}{s}$.
$$
\begin{align}
    \hbar &= \frac{kg m^2}{s} \\
    \hbar &= \frac{kg m^2}{s} \cdot \left( 6.022 \cdot 10^{26} \frac{\mbox{AMU}}{kg} \right) \\
    \hbar &= \frac{\mbox{AMU} m^2}{s} \cdot \left(10^{20} \frac{\mathring A^2}{m^2} \right)\\
    \hbar &= \frac{\mbox{AMU} \mathring A^2}{s} \\
    \\
    \alpha &= \frac{\omega\mu}{\hbar} \\
    \alpha &= \frac{\frac{1}{s}\cdot \mbox{AMU}}{ \left(\frac{\mbox{AMU} \mathring A^2}{s} \right)} \\
    \alpha &= \frac{1}{s}\cdot \mbox{AMU} \cdot \left(\frac{s}{\mbox{AMU} \mathring A^2} \right) \\
    \alpha &= \frac{\mbox{AMU}}{s} \cdot \frac{s}{\mbox{AMU} \mathring A^2} \\
    \alpha &= \frac{1}{\mathring A^2} 
\end{align}
$$
<br>
It should be noted that the central purpose of $\alpha$ is to act as a converter value that makes other values unitless, as well as being used to give the wavefunction final units of $\frac{1}{\sqrt{\mathring A}}$.
<br>
$$
\begin{align}
    y &= \sqrt{\alpha}x \\
    y &= \sqrt{\frac{1}{\mathring A^2}} \cdot \mathring A \\
    y &= \frac{1}{\mathring A} \cdot \mathring A \\
    y &= \frac{\mathring A}{\mathring A} \\
    y &= 1 &&y\mbox{ becomes a unitless quantity.} 
\end{align}
$$
<br>
Below follows a unit analysis of the wavefunction itself taking into consideration the above variable units: 
$$
\begin{align}
    \psi_n &= H_n(y) \cdot e^{-\frac{\alpha x^2}{2}} \\
    \psi_n &= H_n(\mbox{Unitless Quantity}) \cdot e^{-\frac{\frac{1}{\mathring A^2}}{\mathring A^2}} \\
    \psi_n &= \left(\mbox{Unitless Quantity}\right) \cdot e^{-\frac{\mathring A^2}{\mathring A^2}} \\
    \psi_n &= \left(\mbox{Unitless Quantity}\right) \cdot e^{\left(\mbox{Unitless Quantity}\right)} \\
    \psi_n &= \left(\mbox{Unitless Quantity}\right) \cdot \left(\mbox{Unitless Quantity}\right) \\
    \psi_n &= \mbox{Unitless Quantity} 
\end{align}
$$
<br>
As can be seen above, the wavefunction now produces a unitless quantity and all that is left to do is multiply in the unit of $\frac{1}{\sqrt{\mathring A}}$. This is done through $c_0$.
$$
\begin{align}
    c_0 &= \left[\frac{\alpha}{\pi}\right]^\frac{1}{4} \\
    c_0 &= \left[\frac{\frac{1}{\mathring A^2}}{\pi}\right]^\frac{1}{4} \\
    c_0 &= \left[\frac{1}{\mathring A^2}\right]^{\frac{1}{4}} \\
    c_0 &= \frac{1}{\mathring A^{\frac{2}{4}}} \\
    c_0 &= \frac{1}{\mathring A^{\frac{1}{2}}} \\
    c_0 &= \frac{1}{\sqrt{\mathring A}} \\
\end{align}
$$
<br>
Thus, $c_0$ has units of $\frac{1}{\sqrt{\mathring A}}$, meaning that $c_0$ is multipled into the wavefunction to allow the unitless wavefunction to have the correct unit, resulting in the final harmonic oscillator wavefunction: 
$$\psi_n = c_0 \cdot H_n(y) \cdot e^{-\frac{\alpha x^2}{2}}$$


```python
#Allow Notebook to Import from Comp_Chem_Package
import sys
if("win" in sys.platform):
    sys.path.append("..\\..\\Backend")
else: 
    sys.path.append("../../Backend")

from compChemGlobal import *

#Enter Input Data

#Include link to NIST and ask students to find values for a certain diatomic

#Re in units of Angstroms/
#Optimal Diatomic Bond Distance
Re = 1.128323

#w in units of 1/cm
#First Vibrational Constant
w = 2169.81358

#u in units of A
#Reduced Atomic Mass for Diatomic Molecule
u = 12*16/(12+16)

#Anticipate ahead of time? 
#How would expect certain things to change the wavefunctions
#W, Re, and u
#Predict changes in input parameters and 
#appearances of wavefunctions
#Reverse order of notebook, theory first code second
#Very clear for what I am asking
#Scaffold, start with parameters 
#Clear in expectations
#Outcomes and predictiosn
#How developing this code, and all the pourposes 
#Calculate Franc-Condon, transitions
#All in q
```


```python
#All conversion factors acquired from Google

#convert w from 1/cm to 1/m
w *= 100

#convert from 1/m to 1/s using the speed of light, m/s
#multiply by 2pi to adjust waveform
w *= c * 2 * pi

#convert hbar to units of AMU Angstrom^2 / s from kg m^2 / s
#convert from kg to AMU 
hbar2 = hbar
hbar2 *= 6.022 * pow(10, 26)
#convert from m^2 to angstrom^2
hbar2 *= pow(10, 20)

#compute alpha with units of 1/Angstrom^2
alpha = (u*w) / hbar2
```


```python
#Define all functions here
#########################################################################

#converts r to x, where x = r-Re
#r in Angstroms
#Re in Angstroms
#Returns in units of Angstroms
def rToX(r):
    return r - Re

#########################################################################

#Represents the e^x term in the harmonic oscillator
#r in units of angstrom
#unitless return
def eTerm(r):
    x = rToX(r)
    return exp( -alpha * pow(x,2) / 2 )

#########################################################################

#prepares input to the hermite function
#r in units of angstroms
#unitless return
def hermiteInput(r):
    x = rToX(r)
    return sqrt(alpha) * x

#########################################################################

#Function that creates a new Harmonic Oscillator Wavefunction (HOW)
#Each harmonic oscillator function takes a position r in angstroms as input, and returns
#a value in units of 1/sqrt(angstroms)
def newHOW(n):

    normalization = 1 / superSqrt( pow(2,n) * factorial(n) )
    return lambda r : normalization * pow(alpha/pi, 0.25) * eTerm(r) * hermitePolynomials(
        n, hermiteInput(r))
```


```python
#Graphing Code
import numpy as np

fig = plot.go.FigureWidget()
fig.update_layout(
    title_text = "Harmonic Oscillators Wavefunction",
    xaxis_title_text = "r in Angstroms",
    yaxis_title_text = "Wavefunction Output",
)
    
fig.update_layout(showlegend=True)
out = widgets.Output()

basis = []
def graph_wavefunctions(resolution, startR, endR, precision, basisSize, visibleWavefunctions, probability):
    
    visibility = ["legendonly"] * basisSize
    
    try:
        for startEnd in visibleWavefunctions.split(";"):
            if("-" in startEnd):
                startEnd = [int(value) for value in startEnd.split("-")]
                for index in range(startEnd[0], startEnd[1]):
                    visibility[index] = True
            else:
                visibility[int(startEnd)] = True
    except:
        visibility = [False] * basisSize
      
    #build the basis set of the specified size if needed
    if(len(basis) != basisSize):
        for n in range(basisSize):
            basis.append( newHOW(n) )
    
    dr = 1/resolution
    precision = "0." + str(precision)
    
    r = [startR + (dr*step) for step in range( int(abs( (endR - startR)/dr ))) ]
    y = []
    
    name = []
    yCalculations = []
    if("Standard" in probability):
        name.append("HOW")
        yCalculations.append(lambda basisFunction, rValue: basisFunction(rValue)) 
    if("Probability" in probability):
        name.append("HOW Squared")
        yCalculations.append(lambda basisFunction, rValue: basisFunction(rValue) ** 2)
    
    fig.data = []
    for index, basisFunction in enumerate(basis):
        y.append([])
        
        for yCalculation in yCalculations:
            y[index].append([])
            for rValue in r:
                y[index][-1].append( yCalculation(basisFunction, rValue) )
                
        for yIndex, yData in enumerate(y[index]):
            fig.add_trace(plot.go.Scatter(
                name = name[yIndex] + " " + str(index),
                mode = "lines",
                line = {"shape":"spline"},
                x = r, y = yData,

                hovertemplate = 
                    "<b>r = %{x:" + precision + "f}</b><br>" + 
                    "<b>Ψ(r) = %{y:" + precision + "f}</b>",
                hoverlabel = dict(
                    font_size = 16
                ),

                visible = visibility[index]
            ))    
```


```python
#Widgets Code

#code to determine the appropriate start and end values
start = round(-sqrt(-2 * log(pow(10, -15)) / alpha) + Re, 2)
end = round(sqrt(-2 * log(pow(10, -15)) / alpha) + Re, 2)

resolution = widgets.BoundedFloatText(
    value = 100,
    min = .001, 
    step = .01,
    max = pow(10, 10),
    description = '<p style="font-family:verdana;font-size:15px">Resolution</p>',
)

precision = widgets.BoundedIntText(
    value = 2,
    min = 0,
    max = 20,
    step = 1,
    description = '<p style="font-family:verdana;font-size:15px">Precision</p>'
)

basisSize = widgets.BoundedIntText(
    value = 5,
    min = 1, 
    max = pow(10,10),
    step = 1, 
    description = '<p style="font-family:verdana;font-size:15px">Basis Size</p>'
)

visibleWavefunctions = widgets.Text(
    value = "0-" + str(basisSize.value),
    description = '<p style="font-family:verdana;font-size:15px">Visible Ψs</p>'
)

startR = widgets.FloatText(
    value = start, 
    step = .1, 
    description = '<p style="font-family:verdana;font-size:15px">Start r</p>'
)

endR = widgets.FloatText(
    value = end,
    step = .1, 
    description = '<p style="font-family:verdana;font-size:15px">End r</p>'
)

probability = widgets.Dropdown(
    options = ["Standard", "Probability Distribution", "Standard & Probability"],
    description = "Graph Mode"
)

def wavefunctionUpdate(change):
    graph_wavefunctions(resolution.value, startR.value, endR.value, precision.value, basisSize.value, visibleWavefunctions.value, probability.value)
    
resolution.observe(wavefunctionUpdate, "value")
precision.observe(wavefunctionUpdate, "value")
startR.observe(wavefunctionUpdate, "value")
endR.observe(wavefunctionUpdate, "value")
basisSize.observe(wavefunctionUpdate, "value")
visibleWavefunctions.observe(wavefunctionUpdate, "value")
probability.observe(wavefunctionUpdate, "value")

graph_wavefunctions(resolution.value, startR.value, endR.value, precision.value, basisSize.value, visibleWavefunctions.value, probability.value)

display(widgets.VBox([fig, 
             widgets.HBox([resolution, precision, basisSize]),
             widgets.HBox([startR, endR, visibleWavefunctions, probability])]))
display(out)
```


    VBox(children=(FigureWidget({
        'data': [{'hoverlabel': {'font': {'size': 16}},
                  'hovertemplate…



    Output()



```python
import sys
sys.path.append("..\\Comp_Chem_Package")

from how import *
from basisSet import *
from nistScraper import *

dc = getDiatomicConstants("H2")

n1 = how(dc, 20)

bset = basisSet(dc, how, 10)

n1.graph()
bset.graph()
```


    VBox(children=(FigureWidget({
        'data': [{'fill': 'none',
                  'hoverlabel': {'font': {'size': 16}}…



    VBox(children=(FigureWidget({
        'data': [{'fill': 'none',
                  'hoverlabel': {'font': {'size': 16}}…





    <basisSet.basisSet at 0x26e3c6c6b50>




```python
#Building the PES for the Harmonic Oscillator
from howPES import howPES

howP = howPES(dc)
howP.graph()
```


```python
# Solving the Schrodinger Equation for the Harmonic Oscillator 
#Based on equations from the following references: 
#    https://phys.libretexts.org/Bookshelves/University_Physics/Book%3A_University_Physics_(OpenStax)/Map%3A_University_Physics_III_-_Optics_and_Modern_Physics_(OpenStax)/07%3A_Quantum_Mechanics/7.06%3A_The_Quantum_Harmonic_Oscillator#:~:text=Unlike%20a%20classical%20oscillator%2C%20the,%3D%E2%84%8F%CF%89%3Dhf.
#    http://hyperphysics.phy-astr.gsu.edu/hbase/quantum/hosc.html#c1
from operators import *
from basisSet import *
from diatomicConstants import *
from schrodinger import *
from graphable import *
from nistScraper import *

#Declare all variables here
graph = widgets.VBox([])
solution = [0]

#Declare widgets
mass = widgets.FloatText(
    min = pow(10, -8),
    value = 0.5,
    step = .1, 
    description = '<p style="font-family:verdana;font-size:15px">AMU</p>'
)
w = widgets.FloatText(
    min = pow(10, -8),
    value = 2,
    step = 100, 
    description = '<p style="font-family:verdana;font-size:15px">ω in 1/cm</p>'
)
basisSize = widgets.IntText(
    min = 1, 
    value = 5, 
    description = '<p style="font-family:verdana;font-size:15px">Basis Size</p>'
)

dc = buildDiatomicConstants(w=w.value, u=mass.value)
basis = basisSet(dc, size=4)

def buildSolution(graph, returnSolution):
    
    dc["w"] = w.value
    dc["u"] = mass.value
    
    basis = basisSet(dc, size=basisSize.value)
    pes = howPES(dc)

    T = TOperator(basis, integrationStart = -inf)
    V = VOperator(basis, pes, integrationStart = -inf)

    H = T + V
    solution = schrod(H, basis, pes)

    returnSolution.pop()
    returnSolution.append(solution)
    
    #update the graph
    graph.children = solution.graph(getGraph=True).children
    
#Set up widget observation
observeFunction = lambda value : buildSolution(graph, solution)
mass.observe(observeFunction)
w.observe(observeFunction)
basisSize.observe(observeFunction)

#Graph data
buildSolution(graph, solution)
display(graph)
display(widgets.HBox([mass, w, basisSize]))
```


```python
print(solution[0].eigenVectors[0])
```


```python
basis[0]
```


```python

```

fa3aad0c804a8f372b7c66b5f5773c7b0969c5d8
ipynb
Jupyter Notebook
Learning/Basis Sets/Programming HOW.ipynb
GaryZ700/CatLab_CompChem
4a3dee5a3778eb8b69b53b9d0e9ed6321cadedea
MIT
2021-11-13T12:21:52.000Z
2021-11-13T12:21:52.000Z
Learning/Basis Sets/Programming HOW.ipynb
GaryZ700/CatLab_CompChem
4a3dee5a3778eb8b69b53b9d0e9ed6321cadedea
MIT
Learning/Basis Sets/Programming HOW.ipynb
GaryZ700/CatLab_CompChem
4a3dee5a3778eb8b69b53b9d0e9ed6321cadedea
MIT
Qwen/Qwen-72B
1. YES
2. YES