Part-3-Hybrid/Part-3-Hybrid.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4db471a3-a11a-49ab-8318-53821da6a54e",
   "metadata": {},
   "source": [
    "[View on GitHub](https://github.com/csc-training/sscc-cp2k)\n",
    "\n",
    "# Proton transfer simulations with CP2K\n",
    "\n",
    "## Part 3: Re-converging the NEB minimum energy profile using the PBE0 hybrid functional\n",
    "\n",
    "Due to self-interaction error and resulting electron overdelocalization, reaction barriers are typically underestimated by semilocal exchange–correlation functionals such as PBE. Hybrid functionals that include a fraction of exact Hartree–Fock (HF) exchange can mitigate this. We will now re-converge the energies of the NEB images optimized at the PBE level using the **PBE0 hybrid functional**.\n",
    "\n",
    "### Step 3.1: Prepare the input file and submit the calculation\n",
    "\n",
    "**Start by copying the wavefunction restart and coordinate file of each image to the current working directory:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34ab0127-a2af-44ad-885d-9f4a6bb9135e",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cp ../Part-1-NEB/*.wfn .\n",
    "cp ../Part-1-NEB/*Replica* ."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "353e75e0-fcd3-40d9-b0d8-ac3818dd8cfe",
   "metadata": {},
   "source": [
    "CP2K has a built-in `FARMING` program that allows running multiple inputs at the same time within a single job. To re-converge each NEB image, the farming input would look like as follows. **Run the cell to write the input to a file `farming.inp`.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eb5cec0-7399-4de4-b73b-0e82e6268e11",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile farming.inp\n",
    "&GLOBAL\n",
    "  PROJECT farming                         ! Project name, gets prepended to all output files\n",
    "  PROGRAM FARMING                         ! Define farming-type job\n",
    "  RUN_TYPE NONE                           ! Farming master input should have RUN_TYPE NONE\n",
    "&END GLOBAL\n",
    "&FARMING\n",
    "  NGROUPS 8                               ! How many subjobs to run simultaneously\n",
    "  &JOB                                    ! Input files for each subjob follows here\n",
    "    INPUT_FILE_NAME hybrid-1.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-2.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-3.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-4.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-5.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-6.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-7.inp\n",
    "  &END JOB\n",
    "  &JOB\n",
    "    INPUT_FILE_NAME hybrid-8.inp\n",
    "  &END JOB\n",
    "&END FARMING"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7b851d7-e227-430f-8993-b540a7cd043d",
   "metadata": {},
   "source": [
    "As you can see, we need a separate input file for each subjob, in this case one for each image, `hybrid-1.inp` through `hybrid-8.inp`. **Inspect the input file `hybrid-8.inp` and fix all occurrences of `FIXME` based on that:**\n",
    "\n",
    "* We just calculate the energy, so run type should be `ENERGY`\n",
    "* Wavefunction restart file name is `enol-neb-BAND8-RESTART.wfn`\n",
    "* Coordinate file name is `enol-neb-pos-Replica_nr_8-1.xyz`\n",
    "* PBE0 functional *replaces* 25% of the PBE exchange with exact Hartree–Fock (HF) exchange\n",
    "\n",
    "Rest of the inputs are provided and do not need to be edited."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02eb2ac8-19ab-4fbe-841b-67781e2688d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile hybrid-8.inp\n",
    "&GLOBAL\n",
    "  PROJECT enol-hybrid-8                   ! Project name, gets prepended to all output files\n",
    "  RUN_TYPE FIXME                          ! Type of calculation\n",
    "  PRINT_LEVEL LOW                         ! Low verbosity\n",
    "&END GLOBAL\n",
    "\n",
    "&FORCE_EVAL                               ! Include DFT and system settings\n",
    "  &DFT\n",
    "    RESTART_FILE_NAME FIXME\n",
    "    BASIS_SET_FILE_NAME BASIS_MOLOPT\n",
    "    &QS\n",
    "      EPS_DEFAULT 1.0E-12                 ! Sets convergence criteria so that energy should be accurate up to this value\n",
    "    &END QS\n",
    "    &POISSON                              ! How to deal with electrostatics\n",
    "      PERIODIC NONE                       ! System is not periodic\n",
    "      POISSON_SOLVER WAVELET              ! Efficient solver for non-periodic system\n",
    "    &END POISSON\n",
    "    &SCF\n",
    "      SCF_GUESS RESTART                   ! How to guess initial electron density for SCF cycle\n",
    "      EPS_SCF 1.0E-6                      ! Energy convergence criteria (atomic units)\n",
    "      &OT                                 ! Use orbital transformation method instead of brute-force diagonalization\n",
    "        ALGORITHM IRAC                    ! OT algorithm\n",
    "        MINIMIZER DIIS                    ! Minimization algorithm\n",
    "        PRECONDITIONER FULL_ALL           ! Preconditioner for the minimizer\n",
    "      &END OT\n",
    "      MAX_SCF 20                          ! Maximum number of SCF steps before rebuilding preconditioner\n",
    "      &OUTER_SCF ON\n",
    "        EPS_SCF 1.0E-6                    ! Outer loop energy convergence criteria (atomic units)\n",
    "        MAX_SCF 12                        ! Maximum number of preconditioner rebuilds before terminating\n",
    "      &END OUTER_SCF\n",
    "      &PRINT\n",
    "        &RESTART\n",
    "          BACKUP_COPIES 0                 ! Avoid excessive printing of restart files\n",
    "        &END RESTART\n",
    "      &END PRINT\n",
    "    &END SCF\n",
    "    &XC\n",
    "      &XC_FUNCTIONAL\n",
    "        &PBE\n",
    "          SCALE_X FIXME                   ! Include 75% PBE exchange\n",
    "          SCALE_C 1.0                     ! Include 100% PBE correlation\n",
    "        &END PBE\n",
    "      &END XC_FUNCTIONAL\n",
    "      &HF\n",
    "        FRACTION FIXME                    ! Include 25% HF exchange\n",
    "      &END HF\n",
    "    &END XC\n",
    "  &END DFT\n",
    "  &SUBSYS\n",
    "    &CELL\n",
    "      PERIODIC NONE                       ! System is not periodic\n",
    "      ABC 15.0 15.0 15.0                  ! Dimensions of the simulation cell in angstroms\n",
    "    &END CELL\n",
    "    &TOPOLOGY                             ! Specify input coordinate format and filename\n",
    "      COORD_FILE_FORMAT XYZ\n",
    "      COORD_FILE_NAME FIXME\n",
    "    &END TOPOLOGY\n",
    "    &KIND H                               ! Basis sets and pseudopotential for hydrogen\n",
    "      BASIS_SET DZVP-MOLOPT-SR-GTH\n",
    "      POTENTIAL GTH-PBE-q1\n",
    "    &END KIND\n",
    "    &KIND C                               ! Basis sets and pseudopotential for carbon\n",
    "      BASIS_SET DZVP-MOLOPT-SR-GTH\n",
    "      POTENTIAL GTH-PBE-q4\n",
    "    &END KIND\n",
    "    &KIND O                               ! Basis sets and pseudopotential for oxygen\n",
    "      BASIS_SET DZVP-MOLOPT-SR-GTH\n",
    "      POTENTIAL GTH-PBE-q6\n",
    "    &END KIND\n",
    "  &END SUBSYS\n",
    "&END FORCE_EVAL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f5220de-53fd-40bf-a1c7-4659fbcaa930",
   "metadata": {},
   "source": [
    "**Next, write the following simple batch job script to a file by running the cell.** We'll again request two full Puhti nodes (40 CPU cores each), meaning that each of the 8 images will be run using $80/8=10$ cores."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28e98cb8-cda8-4440-9579-477e220f0852",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%writefile job.sh\n",
    "#!/bin/bash\n",
    "#SBATCH --time=00:15:00\n",
    "#SBATCH --reservation=sscc_thu_large\n",
    "#SBATCH --partition=large\n",
    "#SBATCH --nodes=2\n",
    "#SBATCH --ntasks-per-node=40\n",
    "#SBATCH --account=project_2006657\n",
    "\n",
    "module purge\n",
    "module load intel-oneapi-compilers-classic/2021.6.0\n",
    "module load intel-oneapi-mpi/2021.6.0\n",
    "module load cp2k/2024.1\n",
    "\n",
    "srun cp2k.psmp farming.inp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abe8ce5f-29f7-4c3e-a0a8-dc799b7101ec",
   "metadata": {},
   "source": [
    "**Run the calculation as a batch job by executing the cell below.** Note, before running `sbatch`, we unset all Slurm environment variables as we are submitting a batch job from another Slurm job (this notebook that we're currently running)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "059ed3a2-3a16-48f6-933c-7354457deff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "!unset ${!SLURM@}; sbatch job.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68584eb3-57d0-4962-85ba-2bb2daf82912",
   "metadata": {},
   "source": [
    "The calculation should only take about 2 minutes. You may monitor the queue using `squeue` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb1bcec4-9b83-4d4e-ad0a-ec10abeaadd1",
   "metadata": {},
   "outputs": [],
   "source": [
    "!squeue --me"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcf4d265-d5e5-4adf-950d-e7652bc8c2da",
   "metadata": {},
   "source": [
    "### Step 3.2: Analysis\n",
    "\n",
    "**Once the job has completed, extract the PBE0 energies into a file by running the cell below:**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "370af9b9-540f-432d-b40f-74667002f934",
   "metadata": {},
   "outputs": [],
   "source": [
    "!grep \"ENERGY|\" enol-hybrid* | awk '{print $10}' > energy-pbe0.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc86be5e-b3cb-46b8-ac8e-0ca21be94c0e",
   "metadata": {},
   "source": [
    "**Then, plot the PBE0 energy profile together with the PBE profile (in kJ/mol):**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bc14d13-9f5d-472d-a940-384ec2ebdf33",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Read the energy of each image, convert to kJ/mol and shift IS energy to 0\n",
    "energy_pbe = np.loadtxt('../Part-1-NEB/energy-pbe.txt')\n",
    "energy_pbe = 2625.5*(energy_pbe-energy_pbe[0])\n",
    "energy_pbe0 = np.loadtxt('energy-pbe0.txt')\n",
    "energy_pbe0 = 2625.5*(energy_pbe0-energy_pbe0[0])\n",
    "\n",
    "plt.plot(energy_pbe, 'o-', label='PBE')\n",
    "plt.plot(energy_pbe0, 'o-', label='PBE0')\n",
    "plt.xlabel('NEB image')\n",
    "plt.ylabel('Energy (kJ/mol)')\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c3f814a-1728-4fbf-97c8-529bb98bbf39",
   "metadata": {},
   "source": [
    "* Did the barrier height increase compared to the PBE result?\n",
    "\n",
    "**Calculate the equilibrium constant for the proton transfer based on the PBE0 result and:**\n",
    "\n",
    "$$K\\approx\\exp\\left(-\\Delta E/k_\\mathrm{B}T\\right)$$\n",
    "\n",
    "**where $\\Delta E=E_\\mathrm{FS}-E_\\mathrm{IS}$ is the reaction energy.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8596d21a-3c7c-4141-b834-16d0e7d45f66",
   "metadata": {},
   "outputs": [],
   "source": [
    "kT = 2.479  # Boltzmann's constant times 298 K in units of kJ/mol\n",
    "K = np.exp(-energy_pbe0[-1]/kT)\n",
    "print(K)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aac0d1e-c212-4670-8084-6df4d0920942",
   "metadata": {},
   "source": [
    "NMR results suggest that there should be a 76%–24% split between the IS and FS enol forms, respectively. **Estimate this split based on the calculated equilibrium constant.** Assume:\n",
    "\n",
    "$$[\\mathrm{IS}]+[\\mathrm{FS}]=1$$\n",
    "$$\\frac{[\\mathrm{FS}]}{[\\mathrm{IS}]}=K$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c6141f7-ab5c-4c76-83f1-e00b88f4ddf8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve the above system of equations for [IS] and [FS]!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7256f86f-0035-4b7a-a470-cfbfbc951a2b",
   "metadata": {},
   "source": [
    "* How well does your result agree with the experimental observation?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cb82914-5157-40a3-96cf-c353f1b076f3",
   "metadata": {},
   "source": [
    "**[Finish the hands-on with Part 4.](../Part-4-AIMD/Part-4-AIMD.ipynb)**"
   ]
  }
 ],
 "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}