{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "17ff72f2",
   "metadata": {},
   "source": [
    "# 1.2.3 Post-processing Qbox Outputs to compute dynamical matrix.\n",
    "Here we would learn how to post-process the qbox output(s) enmfd-1.r after all the job(s) regarding the displaced coordinates have finished. \n",
    "\n",
    "For that purpose we need to import the following objects from pyepfd."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3329d2d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                 \n",
      "          ███████████                            \n",
      "         ░░███░░░░░███                           \n",
      "          ░███    ░███ █████ ████                \n",
      "          ░██████████ ░░███ ░███                 \n",
      "          ░███░░░░░░   ░███ ░███                 \n",
      "          ░███         ░███ ░███                 \n",
      "          █████        ░░███████                 \n",
      "         ░░░░░          ░░░░░███                 \n",
      "                        ███ ░███                 \n",
      "                       ░░██████                  \n",
      "                        ░░░░░░                   \n",
      " ██████████ ███████████  ███████████ ██████████  \n",
      "░░███░░░░░█░░███░░░░░███░░███░░░░░░█░░███░░░░███ \n",
      " ░███  █ ░  ░███    ░███ ░███   █ ░  ░███   ░░███\n",
      " ░██████    ░██████████  ░███████    ░███    ░███\n",
      " ░███░░█    ░███░░░░░░   ░███░░░█    ░███    ░███\n",
      " ░███ ░   █ ░███         ░███  ░     ░███    ███ \n",
      " ██████████ █████        █████       ██████████  \n",
      "░░░░░░░░░░ ░░░░░        ░░░░░       ░░░░░░░░░░   \n",
      "PyEPFD version     :  1.1\n",
      "Author             : Arpan Kundu\n",
      "Author Email       : arpan.kundu@gmail.com\n",
      "Today              :  2024-12-18 17:19:46.623904\n",
      "*************************************************\n",
      "                   CITATIONS                     \n",
      "=================================================\n",
      "Please cite the following 3 references:          \n",
      "(1) A. Kundu et al, Phys. Rev. Mater (2021), 5,  \n",
      "L070801, \n",
      "(2) A. Kundu and G Galli, \n",
      "J. Chem. Theory. Comput. (2023), 19, 4011\n",
      "(3) https://pyepfd.readthedocs.io/en/latest/\n",
      "*************************************************\n"
     ]
    }
   ],
   "source": [
    "from pyepfd.coord_util import *\n",
    "from pyepfd.elph_classes import *\n",
    "from pyepfd.pyepfd_io import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaa4f350",
   "metadata": {},
   "source": [
    "Like example 2.1, we again read all information from the checkpoint/restart XML file created in the Cartesian phonon calculation. We would instantiate an object named enmfd_inp of class read_pyepfd_info that belongs to coord_util object. The class read_pyepfd_info has several objects such as acoustic sum rule (asr), refined dynamical matrix (ref_dynmatrix), optimized coordinates(coord), atoms and cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2ebb86c8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Process-id0: Time spent on read_pyepfd_info class: 0.0004851818084716797 s.\n"
     ]
    }
   ],
   "source": [
    "enmfd_inp = read_pyepfd_info(\\\n",
    "            file_path='../1_cartesian_phonon/fdphonon.xml')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60aea5cc",
   "metadata": {},
   "source": [
    "Now we would save these informations into different variable, similar to example 2.1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "caee4bc7",
   "metadata": {},
   "outputs": [],
   "source": [
    "asr = enmfd_inp.asr\n",
    "inp_dynmat = enmfd_inp.ref_dynmatrix\n",
    "opt_coord = enmfd_inp.coord\n",
    "atoms = enmfd_inp.atoms\n",
    "cell = enmfd_inp.cell\n",
    "cell_v = abc2h(cell)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "498ef3a6",
   "metadata": {},
   "source": [
    "Now, we would read the qbox output(s) and store them within a list object. Let us first initialize the list object named qbouts. Then we would read the qbox outputs using the qbox class available within coord_utils."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7138a64e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reording atoms, forces, and coordinates.\n",
      "Time spent on qbox class: 0.024749040603637695 s.\n"
     ]
    }
   ],
   "source": [
    "qbouts = []\n",
    "for i in range(1):\n",
    "    qbouts.append(qbox(file_path = 'enmfd-'+str(i+1)+'.r',io = 'r'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fde78a32",
   "metadata": {},
   "source": [
    "Note: during qbox calculation sometimes it rearranges the atom order. Therefore PyEPFD again reverting back to input order. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61b77cae",
   "metadata": {},
   "source": [
    "First, we would read the qbox outputs and store them within a list object. Let us first initialize the list object named qbouts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "13e43714",
   "metadata": {},
   "outputs": [],
   "source": [
    "coords = qbouts[0].coords\n",
    "atoms = qbouts[0].atoms\n",
    "mass = qbouts[0].mass\n",
    "forces = qbouts[0].forces"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11e0cc3a",
   "metadata": {},
   "source": [
    "Next we will append other elements upto the length of the qbouts.(Note as here we have only one qbox output file the next step is unnecessary, however usually we will have more than one qbox output to harness the parallel submission of of multiple jobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "57951ccd",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1,len(qbouts)):\n",
    "    coords = np.vstack((coords,qbouts[i].coords))\n",
    "    forces = np.vstack((forces,qbouts[i].forces))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b57772c1",
   "metadata": {},
   "source": [
    "Now we compute the force constant matrix and diagonaliz it to obtained normal-modes and its frequencies. This is done using the phonon_calculator class available within elph_classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "326e5c5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "phonon = phonon_calculator(forces = forces,\\\n",
    "                           mass = mass,\\\n",
    "                           dynmat = inp_dynmat,\\\n",
    "                           mode = 'ENMFD',\\\n",
    "                           deltax = 0.005,\\\n",
    "                           deltae = 0.001,\\\n",
    "                           ngrid = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "338d7f49",
   "metadata": {},
   "source": [
    "Please remember to use the same mode, deltax and deltae that was used to prepare the input files, i.e., on tutorial enmfd_phonon_1.\n",
    "\n",
    "Next, we would take the computed dynamical matrix and instantiate a dynamica matrix (dm) class available in elph_classes. This class has the methods to refine the dynamical matrix using the suitable acoustic sum rule (asr). Here we have a linear molecule so we will choose asr = 'lin'. For a non-linear polyatomic molecule and crystal, one should use asr = 'lin' & asr = 'crystal', respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "15af4254",
   "metadata": {},
   "outputs": [],
   "source": [
    "dm = dm(dynmat = phonon.dynmat, mass = mass)\n",
    "dm.apply_asr(opt_coord = opt_coord, asr = 'lin')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5debb252",
   "metadata": {},
   "source": [
    "Now it is time to write all information into a checkpoint/restart file that can be used for future reference. This file would be the starting point for calculating a normal-mode finite difference method, stochastic method, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e2222b67",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time spent at write_info       0.0005 s\n"
     ]
    }
   ],
   "source": [
    "restart = write_pyepfd_info(inp_dynmat = inp_dynmat,\\\n",
    "                  dynmat = dm.dynmatrix,\\\n",
    "                  ref_dynmat = dm.refdynmatrix,\\\n",
    "                  mass = mass,\\\n",
    "                  atoms = atoms,\\\n",
    "                  opt_coord = opt_coord,\\\n",
    "                  cell = cell,\\\n",
    "                  file_name='enmfdphonon.xml',\\\n",
    "                  mode='enmfd',\\\n",
    "                  deltax=0.005,\\\n",
    "                  deltae= 0.001,\\\n",
    "                  asr = asr)\n",
    "# Deleting the restart object to finishing the file writing\n",
    "del restart"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c7062d3e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<pyepfd>\n",
      "  <phonon mode='enmfd'>\n",
      "    <ngrid> 1 </ngrid>\n",
      "    <deltax> 0.005 </deltax>\n",
      "    <deltae> 0.001 </deltae>\n",
      "    <asr> lin </asr>\n",
      "    <inp_dynmat shape='(9, 9)'>\n",
      "[  1.16570877e-06, -1.27450455e-12,  1.77449774e-12,  1.16571357e-06, -1.27450979e-12, \n",
      "  -7.58294525e-12, -2.69084475e-06,  2.94198171e-12,  6.70390356e-12, -1.27450455e-12, \n",
      "   1.16571224e-06,  1.48889657e-12, -1.27450979e-12,  1.16571704e-06, -8.50449734e-13, \n",
      "   2.94198171e-12, -2.69085276e-06, -7.36872634e-13,  1.77449774e-12,  1.48889657e-12, \n",
      "   3.44489401e-05,  1.77450504e-12,  1.48890270e-12, -3.16584535e-06, -4.09613279e-12, \n",
      "  -3.43687002e-12, -3.61058355e-05,  1.16571357e-06, -1.27450979e-12,  1.77450504e-12, \n",
      "   1.16571837e-06, -1.27451504e-12, -7.58297646e-12, -2.69085583e-06,  2.94199381e-12, \n",
      "   6.70393115e-12, -1.27450979e-12,  1.16571704e-06,  1.48890270e-12, -1.27451504e-12, \n",
      "   1.16572184e-06, -8.50453234e-13,  2.94199381e-12, -2.69086383e-06, -7.36875666e-13, \n",
      "  -7.58294525e-12, -8.50449734e-13, -3.16584535e-06, -7.58297646e-12, -8.50453234e-13, \n",
      "   3.44501870e-05,  1.75039675e-11,  1.96312172e-12, -3.61072746e-05, -2.69084475e-06, \n",
      "   2.94198171e-12, -4.09613279e-12, -2.69085583e-06,  2.94199381e-12,  1.75039675e-11, \n",
      "   6.21136742e-06, -6.79107529e-12, -1.54748460e-11,  2.94198171e-12, -2.69085276e-06, \n",
      "  -3.43687002e-12,  2.94199381e-12, -2.69086383e-06,  1.96312172e-12, -6.79107529e-12, \n",
      "   6.21138590e-06,  1.70094788e-12,  6.70390356e-12, -7.36872634e-13, -3.61058355e-05, \n",
      "   6.70393115e-12, -7.36875666e-13, -3.61072746e-05, -1.54748460e-11,  1.70094788e-12, \n",
      "   8.33458035e-05 ]\n",
      "    </inp_dynmat>\n",
      "    <dynmat shape='(9, 9)'>\n",
      "[  1.52149567e-06, -2.91302688e-12,  4.40786142e-14,  1.00370699e-06, -1.38621426e-12, \n",
      "  -4.18653074e-12, -2.98797305e-06, -2.22351831e-09,  3.40829809e-12, -2.91302688e-12, \n",
      "   1.45170088e-06,  3.17891915e-12, -9.32830194e-13,  1.08651713e-06,  8.91803023e-13, \n",
      "   2.23137564e-09, -2.96824894e-06,  5.97604638e-12,  4.40786142e-14,  3.17891915e-12, \n",
      "   3.46788383e-05, -2.04652743e-13,  1.37362634e-11, -2.89580040e-06,  2.83697844e-13, \n",
      "   7.15144431e-12, -3.66635112e-05,  1.00370699e-06, -9.32830194e-13, -2.04652743e-13, \n",
      "   1.51822645e-06,  5.27884977e-13, -4.67908801e-12, -2.99049244e-06, -2.22795549e-09, \n",
      "   2.84124381e-12, -1.38621426e-12,  1.08651713e-06,  1.37362634e-11,  5.27884977e-13, \n",
      "   1.33558110e-06,  8.89315199e-12,  2.22812614e-09, -3.02818958e-06,  2.72793760e-11, \n",
      "  -4.18653074e-12,  8.91803023e-13, -2.89580040e-06, -4.67908801e-12,  8.89315199e-12, \n",
      "   3.46791765e-05,  1.04296605e-11,  8.69163359e-12, -3.66638958e-05, -2.98797305e-06, \n",
      "   2.23137564e-09,  2.83697844e-13, -2.99049244e-06,  2.22812614e-09,  1.04296605e-11, \n",
      "   7.14882744e-06, -2.95473689e-12, -6.77145994e-12, -2.22351831e-09, -2.96824894e-06, \n",
      "   7.15144431e-12, -2.22795549e-09, -3.02818958e-06,  8.69163359e-12, -2.95473689e-12, \n",
      "   7.14956898e-06,  1.59883079e-11,  3.40829809e-12,  5.97604638e-12, -3.66635112e-05, \n",
      "   2.84124381e-12,  2.72793760e-11, -3.66638958e-05, -6.77145994e-12,  1.59883079e-11, \n",
      "   8.46401151e-05 ]\n",
      "    </dynmat>\n",
      "    <ref_dynmat shape='(9, 9)'>\n",
      "[  1.31717371e-06, -7.25798794e-13,  6.22844237e-13,  1.31717914e-06, -7.25801781e-13, \n",
      "  -3.76082417e-12, -3.04047638e-06,  1.67538576e-12,  3.62174484e-12, -7.25798794e-13, \n",
      "   1.31717570e-06,  7.02167317e-13, -7.25801781e-13,  1.31718113e-06, -7.55927237e-13, \n",
      "   1.67538576e-12, -3.04048097e-06,  6.20477885e-14,  6.22844237e-13,  7.02167317e-13, \n",
      "   3.46719008e-05,  6.22846800e-13,  7.02170207e-13, -2.90273971e-06, -1.43773230e-12, \n",
      "  -1.62083643e-12, -3.66668359e-05,  1.31717914e-06, -7.25801781e-13,  6.22846800e-13, \n",
      "   1.31718456e-06, -7.25804769e-13, -3.76083964e-12, -3.04048889e-06,  1.67539266e-12, \n",
      "   3.62175975e-12, -7.25801781e-13,  1.31718113e-06,  7.02170207e-13, -7.25804769e-13, \n",
      "   1.31718655e-06, -7.55930348e-13,  1.67539266e-12, -3.04049348e-06,  6.20480439e-14, \n",
      "  -3.76082417e-12, -7.55927237e-13, -2.90273971e-06, -3.76083964e-12, -7.55930348e-13, \n",
      "   3.46722354e-05,  8.68123688e-12,  1.74493226e-12, -3.66672221e-05, -3.04047638e-06, \n",
      "   1.67538576e-12, -1.43773230e-12, -3.04048889e-06,  1.67539266e-12,  8.68123688e-12, \n",
      "   7.01843385e-06, -3.86734929e-12, -8.36019539e-12,  1.67538576e-12, -3.04048097e-06, \n",
      "  -1.62083643e-12,  1.67539266e-12, -3.04049348e-06,  1.74493226e-12, -3.86734929e-12, \n",
      "   7.01844445e-06, -1.43226996e-13,  3.62174484e-12,  6.20477885e-14, -3.66668359e-05, \n",
      "   3.62175975e-12,  6.20480439e-14, -3.66672221e-05, -8.36019539e-12, -1.43226996e-13, \n",
      "   8.46395618e-05 ]\n",
      "    </ref_dynmat>\n",
      "  </phonon>\n",
      "  <cell shape='(6)'>\n",
      "[  2.00000000e+01,  2.00000000e+01,  2.00000000e+01,  9.00000000e+01,  9.00000000e+01, \n",
      "   9.00000000e+01 ]\n",
      "  </cell>\n",
      "  <mass shape='(9)'>\n",
      "[  2.91651223e+04,  2.91651223e+04,  2.91651223e+04,  2.91651223e+04,  2.91651223e+04, \n",
      "   2.91651223e+04,  2.18941669e+04,  2.18941669e+04,  2.18941669e+04 ]\n",
      "  </mass>\n",
      "  <opt_coord shape='(9)'>\n",
      "[  0.00000000e+00,  0.00000000e+00,  2.18666400e+00,  0.00000000e+00,  0.00000000e+00, \n",
      "  -2.18666100e+00,  0.00000000e+00,  0.00000000e+00, -3.00000000e-06 ]\n",
      "  </opt_coord>\n",
      "  <atoms shape='(3)'>\n",
      "[ O, O, C ]\n",
      "  </atoms>\n",
      "</pyepfd>"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cat enmfdphonon.xml"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8952dad7",
   "metadata": {},
   "source": [
    "We see this is an xml file. It is similar to the restart file we created at the end of example 1.1.3, i.e., the one we used as input for normal mode phonon calculations except here is one additional dynmacila matric element named inp_dynamt that is the initial dynamical matrix used to define the normal modes along which displacements are performed. \n",
    "\n",
    "We can visualize the normal modes as we did in example 1.1.3. The below example shows how to do it using VMD file as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "0cea9140",
   "metadata": {},
   "outputs": [],
   "source": [
    "nmodes = write_nmode(atoms = atoms, \n",
    "                     cell_v = cell_v, \n",
    "                     opt_coord = opt_coord, \n",
    "                     mode_v = dm.refV, # refined normal mode vectors after ASR \n",
    "                     mode_freq = dm.refomega, #refined frequencies \n",
    "                     fmt='nmd', #Options are: 'axsf','nmd','molden'\n",
    "                     file_path='refdynmat')\n",
    "# Deleting the nmodes object to finish printing the information into the file\n",
    "del nmodes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07e92b51",
   "metadata": {},
   "source": [
    "This will create an nmd file named refdynmat.nmd cintaining the normal mode vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "92fbefa7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "title PyEPFD-nmodes \n",
      "names O O C \n",
      "coordinates 0.0 0.0 1.157132840007336 0.0 0.0 -1.1571312524755888 0.0 0.0 -1.587531747e-06 \n",
      "mode 1       -44347.6   -0.000245977    -0.00198492    7.05801e-05    -0.00495603    0.000405116    7.05801e-05    -0.00260101   -0.000789901    7.05801e-05 \n",
      "mode 2        -214174    -0.00109769     -0.0020006    0.000150957    0.000918096    -0.00449894    0.000150957   -8.97939e-05    -0.00324977    0.000150957 \n",
      "mode 3         329802    -0.00223329    -0.00398869   -6.44854e-05      0.0020469     0.00299984   -6.44854e-05   -9.31928e-05   -0.000494414   -6.44854e-05 \n",
      "mode 4         101100    -0.00483253     0.00239793   -4.45383e-05    6.49351e-05    -0.00038747   -4.45383e-05    -0.00238379     0.00100523   -4.45383e-05 \n",
      "mode 5        71131.4   -4.99698e-05    8.27297e-05     0.00352578    9.81602e-05    0.000234485     0.00352578    2.40954e-05    0.000158607     0.00352578 \n",
      "mode 6     0.00146652     0.00205646     0.00067057    -1.0886e-09     0.00205647    0.000670573    1.53405e-09    -0.00547883    -0.00178653   -5.93389e-10 \n",
      "mode 7     0.00146652     0.00067057    -0.00205646    7.61687e-11    0.000670573    -0.00205647    6.06901e-11    -0.00178653     0.00547883   -1.82309e-10 \n",
      "mode 8    0.000743305    6.50042e-10    2.16218e-10     0.00414051    6.50045e-10    2.16219e-10    -0.00414049   -1.73184e-09   -5.76047e-10   -2.34076e-08 \n",
      "mode 9    0.000422301    -2.3297e-10   -3.99146e-12     0.00216302   -2.32971e-10   -3.99148e-12     0.00216305    6.20679e-10     1.0634e-11    -0.00576275 "
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "cat refdynmat.nmd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cf7e4fe",
   "metadata": {},
   "source": [
    "This file can be visualized using VMD. If we choose fmt='axsf' / 'molden' that would create .axsf / .molden file respectively. They can be visualized with XCrysden and Molden, respectively."
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}