{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Calculating Spectroscopic Constants from a Potential Energy Surface\"\"\"\n",
    "\n",
    "__authors__ = \"Ashley Ringer McDonald\"\n",
    "__credits__ = [\"Dominic A. Sirianni\"]\n",
    "__email__   = [\"armcdona@calpoly.edu\"]\n",
    "\n",
    "__copyright__ = \"(c) 2008-2019, The Psi4Education Developers\"\n",
    "__license__   = \"BSD-3-Clause\"\n",
    "__date__      = \"2019-11-18\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating Spectroscopic Constants from a Potential Energy Surface  \n",
    "The spectroscopic constants for a diatomic molecule can be calculated from a series of potential energies computed for different bond separations.  In this lab, you will calculate the spectroscopic constants for two molecules and compare their force constants.  You will also study the effect of using the harmonic approximation in determining the vibrational energy levels.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To use this lab: As you work through the lab, be sure to execute each code block.  If you don't, then when you try to write you own code, it will tell you the previous variable are undefined because the prior code blocks were not executed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we import the electronic structure functions into our notebook.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import psi4\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next contruct the molecule using the PSI4 input structure and define the values of the bond length, R.  Here we build an HF molecule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mol_tmpl = \"\"\"H\n",
    "F 1 **R**\"\"\"\n",
    "rvals = [0.8, 0.85, 0.9, 0.95, 1.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the following `psi4.` functions calls uses the Psi4 program.  The `psi4.geometry` function creates the molecule and the `psi4.energy` function calculates the SCF/cc-pVDZ energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "molecules =[]\n",
    "energies = []\n",
    "for r in rvals:\n",
    "    molecule = psi4.geometry(mol_tmpl.replace(\"**R**\", str(r)))\n",
    "    molecules.append(molecule)\n",
    "for mol in molecules:\n",
    "    energy = psi4.energy(\"SCF/cc-pVDZ\", molecule=mol)\n",
    "    energies.append(energy)\n",
    "print(rvals)\n",
    "print(energies)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This makes a nicer chart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for num, r in enumerate(rvals):\n",
    "    print(r, ':', energies[num])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*What is the general trend in the energy for these bond distances? Offer an explaination as to why this seems reasonable.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Student answer box: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Activity: *Plot the energies vs. R*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write you code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next function uses the bonds distances and the energies to derive the spectroscopic constants for the molecule.  These are output to a python dictionary which is called data in this example.  All of the values are given in cm$^{-1}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = psi4.diatomic.anharmonicity(rvals, energies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To access an element of data in a later computation, you use the syntax: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data['nu']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Activity: *Calculate the energies of the n vibrational levels.  Recall that the vibrational energy levels are given by $$ E_n = \\omega_e \\left( n+ \\frac{1}{2} \\right)  - \\omega_e x_e \\left( n + \\frac{1}{2} \\right)^2   $$ where the second term is the anharmonic correction.  Assuming J=0, calculate the energy for the n=0 to n=3 energy levels, with and without the anharmonic correction.  Remember, you can access the $ \\omega_e $ and $ \\omega_e x_e $ values from the `data` dictionary.  Have your code print your values.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Write your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Activity: *Calculate the difference between the harmonic and the anharmonic energies and print these values.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Is it more important to include anharmonicity corrections for ground state energy calculations or excited state calculations? Explain your answer.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write your answer here:  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Activity: *Calculate the $ \\Delta E $ for n=0 to n=1 and n=1 to n=2.  Compare the results for the harmonic approximation and the anharmonic results.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write you code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Is the spacing of the energy levels equal?  Explain your answer.*  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write your answer here: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Activity: *Calculate the force constant for the HF bond in N/m.  $$ \\omega_e = \\frac{1}{2\\pi c} \\sqrt \\frac{k}{\\mu} $$  Remember that $ \\mu $ is the reduced mass of HF.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Activity: *Construct and plot the PES for CO.  Determine $\\omega_e$, $\\omega_e x_e$, $\\mu$, and k for CO.  Add additional cells to your jupyter notebook as needed.* "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: *Based on your results, is anharmonicity a more significant factor for HF or CO?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write your answer here:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: *Is the force constant greater for HF or CO?  Do your results agree with the typical bond orders expected for HF and CO?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Write your answer here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}