Chapter05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ce2d6e30",
   "metadata": {},
   "source": [
    "# Chapter 5. *Ab initio* Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "520d5d2e",
   "metadata": {},
   "source": [
    "*Ab initio* is a Latin term that translates to \"from the beginning.\" In the context of computational chemistry, ab initio methods refer to quantum mechanical calculations performed using the principles of quantum mechanics without any empirical data. These calculations aim to solve the Schrödinger equation for the system of interest from first principles.\n",
    "\n",
    "Ab initio methods find applications in various areas:\n",
    "\n",
    "- Molecular Structure Determination: Determining equilibrium structures of molecules.\n",
    "- Reaction Pathways: Studying the mechanism of chemical reactions.\n",
    "- Spectroscopy: Calculating molecular properties that can be compared with experimental spectra."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "295a0f8f",
   "metadata": {},
   "source": [
    "## 5.1. Hartree-Fock (HF) Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bada76cf",
   "metadata": {},
   "source": [
    "The Hartree-Fock method is one of the fundamental approaches in ab initio calculations. It involves approximating the many-electron wavefunction of a system as a product of single-electron wavefunctions (orbitals).\n",
    "\n",
    "Main Concepts:\n",
    "\n",
    "- Slater Determinant: Represents the antisymmetry property of electron wavefunctions.\n",
    "- Electron Correlation: HF does not consider electron correlation, which leads to the difference between the HF energy and the exact non-relativistic energy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d6e53ba",
   "metadata": {},
   "source": [
    "## 5.2. Electron Correlation Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be689ba0",
   "metadata": {},
   "source": [
    "Post-Hartree-Fock methods aim to include the effects of electron correlation, which HF omits.\n",
    "\n",
    "Popular Methods:\n",
    "\n",
    "- Configuration Interaction (CI)\n",
    "- Coupled Cluster (CC) Theory\n",
    "- Many-Body Perturbation Theory (MBPT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68576094",
   "metadata": {},
   "source": [
    "## 5.3. Basis Sets in *Ab Initio* Calculations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c1bcbc6",
   "metadata": {},
   "source": [
    "The choice of basis set is crucial in ab initio calculations. There are various basis sets tailored for different purposes:\n",
    "\n",
    "- Minimal Basis Sets: e.g., STO-3G\n",
    "- Split-Valence Basis Sets: e.g., 6-31G, 6-311G\n",
    "- Polarization Functions: e.g., 6-31G*\n",
    "- Diffuse Functions: e.g., 6-31+G"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3531086",
   "metadata": {},
   "source": [
    "## 5.4. Advantages and Limitations"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d263639",
   "metadata": {},
   "source": [
    "**Advantages:**\n",
    "\n",
    "- Accuracy: Since ab initio methods are based on first principles, they can provide very accurate results when applied correctly.\n",
    "- Generality: They can be applied to any molecular system without requiring any empirical data.\n",
    "\n",
    "**Limitations:**\n",
    "\n",
    "- Computational Cost: Ab initio methods can be computationally expensive, especially for larger molecules or systems.\n",
    "- Basis Set Dependence: The accuracy of results often depends on the choice of basis set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba77f956",
   "metadata": {},
   "source": [
    "## 5.5. Self-Consistent Field (SCF) Method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "092ec747",
   "metadata": {},
   "source": [
    "The Self-Consistent Field (SCF) method is a foundational concept in ab initio calculations, particularly in the context of the Hartree-Fock (HF) theory and its extensions. SCF is a computational approach used to approximate the electronic structure of molecules in quantum chemistry.\n",
    "Key Aspects of the SCF Method:\n",
    "\n",
    "- **Iterative Procedure:** The SCF method employs an iterative process to solve the quantum mechanical equations. Starting with an initial guess of the wavefunction, it refines this guess in a series of steps to achieve a consistent solution that satisfies the Schrödinger equation for the system.\n",
    "\n",
    "- **Fock Matrix and Density Matrix:** The method involves constructing and solving the Fock matrix, which depends on the electron density. The electron density is represented by the density matrix, which is updated in each iteration based on the current wavefunction.\n",
    "\n",
    "- **Convergence Criteria:** The iterative process continues until the difference between successive approximations of the wavefunction (or the total energy) falls below a pre-defined threshold, indicating convergence to a self-consistent solution.\n",
    "- **Role in Hartree-Fock Theory:** In Hartree-Fock theory, the SCF method is used to find the optimal set of molecular orbitals that minimize the total electronic energy of the system. This approach considers the mean-field approximation, where each electron moves in an average field created by all other electrons.\n",
    "\n",
    "***Applications and Significance:***\n",
    "\n",
    "- Electronic Structure Analysis: SCF is a critical step in determining the electronic structure of molecules and forms the basis for more complex computational methods, such as post-Hartree-Fock correlation techniques.\n",
    "- Foundational for Advanced Methods: Methods like Density Functional Theory (DFT) and post-Hartree-Fock methods also use the SCF approach as a part of their computational process.\n",
    "\n",
    "***Limitations:***\n",
    "\n",
    "- Mean-Field Approximation: While the SCF method provides a robust framework, it relies on the mean-field approximation, which does not account for electron correlation effects adequately. This is particularly notable in the basic Hartree-Fock implementation.\n",
    "- Initial Guess Dependency: The final converged solution can sometimes depend on the initial guess, particularly in systems with near-degenerate electronic states."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76216bda",
   "metadata": {},
   "source": [
    "## 5.6. Running *Ab initio* Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38a6d324",
   "metadata": {},
   "source": [
    "We will run *Ab initio* calculations using `psi4`, an open-source computation chemistry framework. To install `psi4`, run the following command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72aa2634",
   "metadata": {},
   "outputs": [],
   "source": [
    "!conda install psi4 -c conda-forge/label/libint_dev"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0009d26e",
   "metadata": {},
   "source": [
    "### 5.6.1. Energy Calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a326b8f",
   "metadata": {},
   "source": [
    "In the following section, we will calculate the energy of 1,3-butadiene conformers using Hartree-Fock method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7dc680d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from rdkit import Chem\n",
    "from rdkit.Chem import AllChem\n",
    "from utils import View3DModel\n",
    "import psi4\n",
    "import py3Dmol\n",
    "from tqdm import tqdm # Used to make progress bar, run `!pip install tqdm`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cfdb69a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a molecule of 1,3-butadiene\n",
    "mol = Chem.MolFromSmiles('C=CC=C')\n",
    "\n",
    "# Prepare the molecule\n",
    "mol = Chem.AddHs(mol)\n",
    "Chem.rdDistGeom.EmbedMolecule(mol, useRandomCoords=True)\n",
    "AllChem.UFFOptimizeMolecule(mol, maxIters=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86278bdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# View 3D model with py3Dmol\n",
    "View3DModel(mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e137353",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the geometry to XYZ string\n",
    "xyz_string = Chem.MolToXYZBlock(mol)\n",
    "\n",
    "# Get the psi4 geometry\n",
    "geometry = psi4.geometry(xyz_string)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c9bfc0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set calculation options\n",
    "psi4.set_options({'basis': 'STO-3G'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8257236",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform Hartree-Fock calculation to compute the energy\n",
    "energy = psi4.energy('HF') * psi4.constants.hartree2kcalmol # Convert from Hartrees to kcal/mol\n",
    "\n",
    "# Print the calculated energy\n",
    "print(f\"The Hartree-Fock energy of 1,3-butadiene with the 3-21G basis set is: {energy:.6f} kcal/mol\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8090e41",
   "metadata": {},
   "source": [
    "Next, we will rotate the molecule around the C2-C3 bond, and calculate the energy for each comformers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18d79270",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defind a function that calculate HF energy\n",
    "def calculate_HF_energy(mol, rotation_angle):\n",
    "    # Rotate round C-C bond\n",
    "    Chem.rdMolTransforms.SetDihedralDeg(mol.GetConformer(0), 0, 1, 2, 3, -180 + rotation_angle)\n",
    "    \n",
    "    # Write the geometry to XYZ string\n",
    "    xyz_string = Chem.MolToXYZBlock(mol)\n",
    "\n",
    "    # Get the psi4 geometry\n",
    "    geometry = psi4.geometry(xyz_string)\n",
    "    \n",
    "    # Calculate the energy total energy using Hartree-Fock\n",
    "    energy = psi4.energy('HF') * psi4.constants.hartree2kcalmol # Convert from Hartrees to kcal/mol\n",
    "    \n",
    "    return energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df7ab3f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate an array of energies for all comformers\n",
    "energy_arr = []\n",
    "angles = np.linspace(0, 360, num=121) # rotate 3 degrees each time\n",
    "progress_bar = tqdm(angles)\n",
    "for angle in progress_bar:\n",
    "    energy_arr.append(calculate_HF_energy(mol, angle))\n",
    "    progress_bar.set_description(f'Current angle: {angle} degree')\n",
    "energy_arr = np.array(energy_arr)\n",
    "    \n",
    "# Plot the result\n",
    "plt.plot(angles, energy_arr)\n",
    "plt.xlabel('Dihedral angle (degree)')\n",
    "plt.ylabel('Energy (kcal/mol)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80b7501c",
   "metadata": {},
   "source": [
    "### 5.6.2. Geometry Optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cd81776",
   "metadata": {},
   "source": [
    "In the following section, we will optimize the geometry of a distorted D-glucose molecule using SCF method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da2497ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "mol = Chem.MolFromPDBFile('./structures/glucose.pdb')\n",
    "\n",
    "# Prepare the molecule\n",
    "mol = Chem.AddHs(mol)\n",
    "Chem.rdDistGeom.EmbedMolecule(mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7fc82ca0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# View 3D model with py3Dmol\n",
    "View3DModel(mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d02c1daa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write the geometry to XYZ string\n",
    "xyz_string = Chem.MolToXYZBlock(mol)\n",
    "\n",
    "# Get the psi4 geometry\n",
    "geometry = psi4.geometry(xyz_string)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7e66c8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set calculation options\n",
    "psi4.set_options({'BASIS': 'STO-3G'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "269ac615",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run SCF geometry optimization\n",
    "max_iters = 500\n",
    "_, wfn_history = psi4.optimize('SCF', molecule=geometry, optking__geom_maxiter=max_iters, return_history=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9b09dbb",
   "metadata": {},
   "source": [
    "Now we can view the optimization result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c1b6bd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the energy values\n",
    "energy_values = wfn_history['energy']\n",
    "plt.plot(energy_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9aa932e7",
   "metadata": {},
   "source": [
    "We can also view the conformers during the optimization process, or even make an animation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f32b28fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get the list of coordinates\n",
    "coordinates_list = wfn_history['coordinates']\n",
    "\n",
    "# Get the number of conformers\n",
    "num_confs = len(coordinates_list)\n",
    "print(num_confs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1137bbea",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function to view the 3D model of a conformer\n",
    "def show_geometry(mol, index):\n",
    "    coords = coordinates_list[index]\n",
    "    \n",
    "    # Adjust the RDKit molecule's atom positions to match the current geometry\n",
    "    conf = mol.GetConformer()\n",
    "    for i, coord in enumerate(coords):\n",
    "        conf.SetAtomPosition(i, coord)\n",
    "    \n",
    "    # Show the conformer\n",
    "    View3DModel(mol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d33f5349",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# For example, visualize the last geometry:\n",
    "show_geometry(mol, num_confs - 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d8e5074",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function to generate a PDB trajectory\n",
    "def generate_pdb_trajectory(mol, coordinates_list):\n",
    "    # Initialize the trajectory string\n",
    "    trajectory = \"\"\n",
    "        \n",
    "    # Loop over each set of coordinates and generate the corresponding PDB format\n",
    "    for idx, coordinates in enumerate(coordinates_list):\n",
    "        # Adjust the molecule's atom positions to match the current geometry\n",
    "        conf = mol.GetConformer()\n",
    "        for i, coord in enumerate(coordinates):\n",
    "            conf.SetAtomPosition(i, coord)\n",
    "            \n",
    "        # Add the records for atoms\n",
    "        trajectory += \"MODEL     {:4d}\\n\".format(idx + 1)\n",
    "        trajectory += Chem.MolToPDBBlock(mol)\n",
    "        trajectory += \"ENDMDL\\n\"\n",
    "\n",
    "    return trajectory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce4ab6f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the PDB trajectory string\n",
    "pdb_trajectory = generate_pdb_trajectory(mol, coordinates_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10a71dc2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Show the animation\n",
    "view = py3Dmol.view(width=400, height=300)\n",
    "view.addModelsAsFrames(pdb_trajectory, \"trajectory\", {'keepH': True})\n",
    "view.setBackgroundColor('white')\n",
    "view.setStyle({'stick': {'scale': 0.3}, 'sphere': {'scale': 0.3}})\n",
    "view.zoomTo()\n",
    "view.animate({'loop': \"forward\", 'interval': 200}) # Adjust the speed as needed (set 'interval' to a new value in millisecond)\n",
    "view.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13f5450b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  },
  "toc": {
   "base_numbering": "",
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Table of Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "203px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}