notebooks/Case_study_large_peptide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "564e4861",
   "metadata": {},
   "source": [
    "# Calculate DFT/CPCM-level ion-pair interaction energies of 4 peptides"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7e93b9a",
   "metadata": {},
   "source": [
    "The following code extracts the total energies of the species A, B, C and D (defined in Fig.2 or Ref.1) from the output files of Orca and calculates the ion-pair interaction energies, $E_{\\rm ip}$, according to Eq.3 of Ref.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0e20fec4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import re\n",
    "import numpy as np\n",
    "au2kcal = 627.50957099203276"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b998cfce",
   "metadata": {},
   "outputs": [],
   "source": [
    "pattern='FINAL SINGLE POINT ENERGY'\n",
    "position=4\n",
    "\n",
    "mols=[\n",
    "    '../Case_Study/1EJG/new_frag/',\n",
    "    '../Case_Study/1BDK/new_frag/',\n",
    "    '../Case_Study/1L2Y/new_frag/',\n",
    "    '../Case_Study/1SCO/new_frag/'\n",
    "]\n",
    "\n",
    "Ctypes=np.array([\n",
    "    [3,2],\n",
    "    [3,2],\n",
    "    [3,2],\n",
    "    [3,2]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbb20c50",
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_E(filename):\n",
    "    file = open(filename, \"r\")\n",
    "    iline=0\n",
    "    for line in file:\n",
    "        if re.search(pattern, line):\n",
    "            iline=iline+1\n",
    "            eline=line\n",
    "    if (iline < 1):\n",
    "        print('ERROR: Can\\'t find <',pattern,'> in ',filename)\n",
    "    file.close()\n",
    "    E=float( eline.split()[position] )\n",
    "    return E\n",
    "\n",
    "def calc_E_SB(mol,Ctype):\n",
    "    \n",
    "    dire=mol\n",
    "    \n",
    "    file_A_pos=dire+'A_pos.out'\n",
    "    file_A_neg=dire+'A_neg.out'\n",
    "    \n",
    "    E_A_pos=extract_E(file_A_pos)\n",
    "    E_A_neg=extract_E(file_A_neg)\n",
    "    A=E_A_pos+E_A_neg\n",
    "    \n",
    "    file_B=dire+'B.out'\n",
    "\n",
    "    B=extract_E(file_B)\n",
    "\n",
    "    N_D=np.prod(Ctype)\n",
    "    E=np.zeros(N_D)\n",
    "    \n",
    "    i_D=0\n",
    "    \n",
    "    for i_Cm in range(Ctype[1]):\n",
    "        for i_Cp in range(Ctype[0]):\n",
    "        \n",
    "            file_C_pos=dire+'C'+str(i_Cp+1)+'_pos.out'\n",
    "            file_C_neg=dire+'C'+str(i_Cm+1)+'_neg.out'\n",
    "        \n",
    "            E_C_pos=extract_E(file_C_pos)\n",
    "            E_C_neg=extract_E(file_C_neg)        \n",
    "            C=E_C_pos+E_C_neg\n",
    "        \n",
    "            file_D=dire+'D'+str(i_D+1)+'.out'\n",
    "        \n",
    "            D=extract_E(file_D)\n",
    "        \n",
    "            E[i_D]=((B - D) - (A - C))*au2kcal\n",
    "            \n",
    "            i_D=i_D+1\n",
    "        \n",
    "    return E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "00fd7638",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Individial ion-pair energies of peptide=    0 \n",
      "[-24.97677648 -25.99801072 -26.09987677 -27.13293271 -19.67280995\n",
      " -24.38353481]\n",
      "Individial ion-pair energies of peptide=    1 \n",
      "[-4.75268083 -4.61003908 -6.32199699 -2.55172504 -2.34532343 -7.15575426]\n",
      "Individial ion-pair energies of peptide=    2 \n",
      "[-6.14369516 -2.90393668 -4.15987735 -4.75146263 -2.13520878 -3.09752916]\n",
      "Individial ion-pair energies of peptide=    3 \n",
      "[ 1.63733782  2.35477589  2.41031416 -1.92283748  1.5349265   1.69050076]\n"
     ]
    }
   ],
   "source": [
    "Nmol=len(mols)\n",
    "\n",
    "Eip=np.zeros([Nmol])\n",
    "Eip_fluct=np.zeros([Nmol])\n",
    "  \n",
    "for imol in range(Nmol):\n",
    "    mol=mols[imol]\n",
    "    Ctype=Ctypes[imol]\n",
    "    E=calc_E_SB(mol,Ctype)\n",
    "    fstring=(f'''Individial ion-pair energies of peptide={imol:5d} ''')\n",
    "    print(fstring)\n",
    "    fstring=(E) \n",
    "    print(fstring)\n",
    "    Eip[imol]=np.mean(E)\n",
    "    Eip_fluct[imol]=np.std(E)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bb6bc74a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Peptide= 0   E_IP_Mean=  -24.71   E_IP_StdDeviation=    2.42\n",
      "Peptide= 1   E_IP_Mean=   -4.62   E_IP_StdDeviation=    1.77\n",
      "Peptide= 2   E_IP_Mean=   -3.87   E_IP_StdDeviation=    1.33\n",
      "Peptide= 3   E_IP_Mean=    1.28   E_IP_StdDeviation=    1.47\n"
     ]
    }
   ],
   "source": [
    "for imol in range(Nmol):\n",
    "    fstring=(f'''Peptide={imol:2d}   E_IP_Mean={Eip[imol]:8.2f}   E_IP_StdDeviation={Eip_fluct[imol]:8.2f}''') \n",
    "    print(fstring)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba172597",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}