General HEOM code for sum of exponentials

heom_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# General heirarchy model/method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qutip.models.heirarchy import Heirarchy\n",
    "from qutip import sigmaz, sigmax, basis, Qobj\n",
    "from qutip.solver import Options\n",
    "import numpy as np\n",
    "\n",
    "from math import pi\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exponentials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ckAR=[-0.0531994+0.359753j, -(0.0531994+0.359753j),(0.265019-0.737435j), (0.265019+0.737435j)]\n",
    "\n",
    "vkAR=[-2.02381-3.37591j,-2.02381+3.37591j,-1.90802-1.16867j,-1.90802+1.16867j]\n",
    "vkAR=[-vk for vk in vkAR]    \n",
    "\n",
    "ckAI=[0.0571124-0.0269892j, 0.0571124+0.0269892j,-0.0573738+0.00965229j, -(0.0573738+0.00965229j)]      \n",
    "vkAI=[-1.70538-3.55923j,-1.70538+3.55923j,-1.40539-1.46163j,-1.40539+1.46163j]\n",
    "vkAI=[-vk for vk in vkAI]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# System specifications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "NR=2\n",
    "NI=2\n",
    "\n",
    "Q = sigmaz()\n",
    "Del =np.pi/2.    \n",
    "wq= 0.0     # Energy of the 2-level system.\n",
    "Hsys = 0.5 * wq * sigmaz() + 0.5 * Del * sigmax()\n",
    "\n",
    "wc = 0.05                # Cutoff frequency.\n",
    "alpha = 0.1/pi           # Coupling strength.\n",
    "N = 15                   # Number of cavity fock states.\n",
    "Temperature = 1/0.95     # Tempertaure.\n",
    "tlist = np.linspace(0, 1, 100)\n",
    "\n",
    "#for amirs Omega= pi data, tlist is different (actually a bit weird)\n",
    "#tlist=tlistA\n",
    "\n",
    "initial_state = basis(2,1) * basis(2,1).dag()                # Initial state of the system.\n",
    "#return_vals = [tensor(qeye(N), kk) for kk in [Q]]            # List for which to calculate expectation value\n",
    "return_vals = [initial_state, basis(2,0) * basis(2,1).dag()   ] # List for which to calculate expectation value\n",
    "eigen_sparse = False\n",
    "calc_time = True                                             \n",
    "options = Options(nsteps=15000, store_states=True,rtol=1e-17, atol=1e-17) # Options for the solver.\n",
    "\n",
    "#Convergence parameters\n",
    "#Nc is the cut-off, and depends on Nc>> \n",
    "Nc = 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Invoking the model and the solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Heirarchy(hamiltonian=Hsys, coupling=Q, Nc=Nc,\n",
    "                  real_prefactors= ckAR, real_exponents=vkAR,\n",
    "                  complex_prefactors=ckAI, complex_exponents=vkAI)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "liouv = Qobj(model.liouvillian())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Quantum object: dims = [[660], [660]], shape = (660, 660), type = oper, isherm = False\\begin{equation*}\\left(\\begin{array}{*{11}c}0.0 & -0.785j & 0.785j & 0.0 & 0.0 & \\cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\-0.785j & 0.0 & 0.0 & 0.785j & 0.0 & \\cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\0.785j & 0.0 & 0.0 & -0.785j & 0.0 & \\cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.785j & -0.785j & 0.0 & 0.0 & \\cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\(-0.476-0.080j) & 0.0 & 0.0 & 0.0 & (-1.405+1.462j) & \\cdots & 0.0 & 0.0 & 0.0 & 0.0 & 0.0\\\\\\vdots & \\vdots & \\vdots & \\vdots & \\vdots & \\ddots & \\vdots & \\vdots & \\vdots & \\vdots & \\vdots\\\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \\cdots & (-6.071-3.376j) & 0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \\cdots & 0.0 & (-6.071-10.128j) & -0.785j & 0.785j & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \\cdots & 0.0 & -0.785j & (-6.071-10.128j) & 0.0 & 0.785j\\\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \\cdots & 0.0 & 0.785j & 0.0 & (-6.071-10.128j) & -0.785j\\\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & \\cdots & 0.0 & 0.0 & 0.785j & -0.785j & (-6.071-10.128j)\\\\\\end{array}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "Quantum object: dims = [[660], [660]], shape = (660, 660), type = oper, isherm = False\n",
       "Qobj data =\n",
       "[[ 0.00000 +0.j          0.00000 -0.78539816j  0.00000 +0.78539816j ...,\n",
       "   0.00000 +0.j          0.00000 +0.j          0.00000 +0.j        ]\n",
       " [ 0.00000 -0.78539816j  0.00000 +0.j          0.00000 +0.j         ...,\n",
       "   0.00000 +0.j          0.00000 +0.j          0.00000 +0.j        ]\n",
       " [ 0.00000 +0.78539816j  0.00000 +0.j          0.00000 +0.j         ...,\n",
       "   0.00000 +0.j          0.00000 +0.j          0.00000 +0.j        ]\n",
       " ..., \n",
       " [ 0.00000 +0.j          0.00000 +0.j          0.00000 +0.j         ...,\n",
       "  -6.07143-10.12773j     0.00000 +0.j          0.00000 +0.78539816j]\n",
       " [ 0.00000 +0.j          0.00000 +0.j          0.00000 +0.j         ...,\n",
       "   0.00000 +0.j         -6.07143-10.12773j     0.00000 -0.78539816j]\n",
       " [ 0.00000 +0.j          0.00000 +0.j          0.00000 +0.j         ...,\n",
       "   0.00000 +0.78539816j  0.00000 -0.78539816j -6.07143-10.12773j   ]]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "liouv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = model.solve(initial_state, tlist, options=options)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}