{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating NQR frequencies\n",
    "===========================\n",
    "\n",
    "In this, we will demonstrate how to calculate NQR frequencies from an ab-initio EFG calculation. This is based on the nqr.py script, which implements a command-line utility to perform the below.\n",
    "\n",
    "You can use this utility from the command line like:\n",
    "\n",
    "```bash\n",
    "> nqr.py [species] [directory]\n",
    "```\n",
    "\n",
    "e.g. if you're in a calculation directory and you want the chlorine couplings you do\n",
    "\n",
    "```bash\n",
    "> nqr.py Cl .\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definition\n",
    "----------\n",
    "\n",
    "From \"Nuclear Quadrupole Resonance Spectroscopy\" by Hand and Das, Solid State Physics Supplement 1, the NQR transition frequency $\\omega$ for a spin $S$ with electric field gradient principal component $V_{zz}$ and quadrupole moment $Q$ is\n",
    "\n",
    "\n",
    "$\\begin{align}\n",
    "    A &= \\frac{eV_{zz}Q}{4S(2S-1)}\\\\\n",
    "    \\omega &= 3 A (2 |m| + 1)/\\hbar\n",
    "\\end{align}$\n",
    "\n",
    "for all available transitions $m \\rightarrow m+1$ where $-S \\le m \\le S$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The calculation\n",
    "---------------\n",
    "\n",
    "We do an EFG calculation using CASTEP with the following input files:\n",
    "\n",
    "### NaClO3.param\n",
    "\n",
    "```raw\n",
    "xcfunctional: PBE \n",
    "opt_strategy: speed\n",
    "task: magres\n",
    "magres_task: efg\n",
    "cut_off_energy: 50 ry\n",
    "```\n",
    "\n",
    "### NaClO3.cell\n",
    "\n",
    "```raw\n",
    "%block SPECIES_POT\n",
    "Pb  3|2.4|2.35|1.6|9.2|12.9|16.5|60UU:61UU:52UU[]\n",
    "%endblock SPECIES_POT\n",
    "\n",
    "%block POSITIONS_ABS\n",
    "ang\n",
    "O 2.664659 0.007779 1.298116\n",
    "O 3.295805 1.990002 3.910951\n",
    "O 0.007779 1.298116 2.664659\n",
    "O 1.298116 2.664659 0.007779\n",
    "O 4.585891 0.622991 6.567871\n",
    "O 1.990002 3.910951 3.295805\n",
    "O 3.910944 3.295798 1.990008\n",
    "O 3.280141 5.277909 5.952599\n",
    "O 0.622985 6.567871 4.585885\n",
    "O 5.277909 5.952599 3.280141\n",
    "O 6.567871 4.585885 0.622991\n",
    "O 5.952599 3.280141 5.277909\n",
    "Na 0.445438 0.445438 0.445438\n",
    "Na 3.733280 2.842423 6.130271\n",
    "Na 2.842423 6.130271 3.733280\n",
    "Na 6.130271 3.733280 2.842423\n",
    "Cl 3.833830 6.030029 0.545995\n",
    "Cl 2.742024 2.742024 2.742024\n",
    "Cl 0.545995 3.833823 6.030023\n",
    "Cl 6.030029 0.545995 3.833830\n",
    "%endblock POSITIONS_ABS\n",
    "\n",
    "%block LATTICE_CART\n",
    "ang\n",
    "6.575801 0.000000 0.000000\n",
    "0.000000 6.575801 0.000000\n",
    "0.000000 0.000000 6.575801\n",
    "%endblock LATTICE_CART\n",
    "\n",
    "kpoint_mp_grid 1 1 1\n",
    "kpoint_mp_offset 0.25 0.25 0.25\n",
    "```\n",
    "\n",
    "and, if just running locally rather than using a cluster submission system, run like\n",
    "\n",
    "```bash\n",
    "> castep.serial NaClO3\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Code\n",
    "----\n",
    "\n",
    "First, we need to make a few imports to get the magres atoms datastructure, constants for units conversions and maths libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "from magres.atoms import MagresAtoms\n",
    "from magres.constants  import millibarn, megahertz\n",
    "import numpy, math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's define a function that will tell us the starting *m* for all available transitions *m->m+1* for a particular nuclear spin *S*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S=0.5, available |m| =  []\n",
      "S=1.0, available |m| =  [0.0]\n",
      "S=1.5, available |m| =  [0.5]\n",
      "S=2.0, available |m| =  [0.0, 1.0]\n"
     ]
    }
   ],
   "source": [
    "def available_m(spin):\n",
    "     return [m for m in numpy.arange(-spin, spin+1, 1) if m >= 0.0][:-1]\n",
    "\n",
    "print (\"S=0.5, available |m| = \", available_m(0.5))\n",
    "print (\"S=1.0, available |m| = \", available_m(1.0))\n",
    "print (\"S=1.5, available |m| = \", available_m(1.5))\n",
    "print (\"S=2.0, available |m| = \", available_m(2.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And let's define a function that will tell us the NQR transition frequency for a particular atom and starting *m*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def calc_nqr(atom, m):\n",
    "    efg = atom.efg\n",
    "    Q = atom.Q\n",
    "    spin = atom.spin\n",
    "    \n",
    "    Vzz = efg.evals[2]\n",
    "    vec_zz = efg.evecs[2]\n",
    "    eta = (abs(efg.evals[0]) - abs(efg.evals[1]))/efg.evals[2]\n",
    "    A = Vzz * (Q * millibarn) / (4.0 * spin * (2.0*spin - 1.0))\n",
    "    fq = 3*A * (2.0*abs(m) + 1.0) * math.sqrt(1.0 + eta**2/3)\n",
    "\n",
    "    return fq\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these we can now loop over all chlorine nuclei in our system and calculate their NQR frequencies for all available transitions. Let's import our ab-initio EFG calculation's output .magres file,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "atoms = MagresAtoms.load_magres('../samples/NaClO3.magres')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and then perform the calculation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35Cl1 S=1.5 Q=-81.65 millibarn\n",
      "  m=0.5->1.5 freq = 28.844 MHz\n",
      "35Cl2 S=1.5 Q=-81.65 millibarn\n",
      "  m=0.5->1.5 freq = 29.057 MHz\n",
      "35Cl3 S=1.5 Q=-81.65 millibarn\n",
      "  m=0.5->1.5 freq = 28.844 MHz\n",
      "35Cl4 S=1.5 Q=-81.65 millibarn\n",
      "  m=0.5->1.5 freq = 28.844 MHz\n",
      "Mean freq = 28.898 Mhz\n"
     ]
    }
   ],
   "source": [
    "freqs = []\n",
    "\n",
    "for Cl_atom in atoms.species('Cl'):\n",
    "    print(Cl_atom, \"S={:.1f} Q={:.2f} millibarn\".format(Cl_atom.spin, Cl_atom.Q))\n",
    "    \n",
    "    for m in available_m(Cl_atom.spin):\n",
    "        freq_MHz = calc_nqr(Cl_atom, m) / megahertz\n",
    "        print(\"  m={:.1f}->{:.1f} freq = {:.3f} MHz\".format(m, m+1, freq_MHz))\n",
    "        \n",
    "        freqs.append(freq_MHz)\n",
    "        \n",
    "print(\"Mean freq = {:.3f} Mhz\".format(numpy.mean(freqs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean frequency is close to the experimental value of ~30 Mhz."
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}