radicalpair-1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import functools\n",
    "import numpy as np\n",
    "import scipy as sci\n",
    "import matplotlib.pyplot as plt\n",
    "import qutip as qt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This dictionary maps string keys ('x', 'y', 'z', 'p', 'm', 'i') to functions that generate spin operators for a given dimension dim.\n",
    "opstr2fun = {'x': lambda dim: qt.spin_Jx((dim-1)/2),\n",
    "             'y': lambda dim: qt.spin_Jy((dim-1)/2),\n",
    "             'z': lambda dim: qt.spin_Jz((dim-1)/2),\n",
    "             'p': lambda dim: qt.spin_Jp((dim-1)/2),\n",
    "             'm': lambda dim: qt.spin_Jm((dim-1)/2),\n",
    "             'i': qt.identity}\n",
    "# Initializes ops as a list of identity matrices for each dimension in dims. Iterates over specs to replace the identity matrix at the specified index with the corresponding spin operator. Returns the tensor product of the operators in ops using qt.tensor.\n",
    "def mkSpinOp(dims, specs):\n",
    "    ops = [qt.identity(d) for d in dims]\n",
    "    for ind, opstr in specs:\n",
    "        ops[ind] = ops[ind] * opstr2fun[opstr](dims[ind])\n",
    "    return qt.tensor(ops)\n",
    "# Constructs a Hamiltonian for a single spin system with interactions along the x, y, and z axes.\n",
    "def mkH1(dims, ind, parvec):\n",
    "    axes = ['x', 'y', 'z']\n",
    "    # Creates a list of spin operators weighted by the corresponding parameters in parvec (ignores zero parameters). Uses functools.reduce to sum these weighted spin operators.\n",
    "    return functools.reduce(lambda a, b: a + b, \n",
    "               [v * mkSpinOp(dims, [(ind,ax)]) for v, ax in zip(parvec, axes) if v!=0])\n",
    "# Constructs a Hamiltonian for the interaction between two spin systems with interaction terms along all combinations of x, y, and z axes.\n",
    "def mkH12(dims, ind1, ind2, parmat):\n",
    "    axes = ['x', 'y', 'z']\n",
    "    ops = []\n",
    "    # Iterates over all combinations of the x, y, and z axes for the two spins. For each non-zero element in parmat, adds the corresponding spin-spin interaction term to the empty list ops.\n",
    "    for i in range(3):\n",
    "        for j in range(3):\n",
    "            if parmat[i,j] != 0:\n",
    "                ops.append(parmat[i,j] * mkSpinOp(dims, [(ind1,axes[i]), (ind2,axes[j])]))\n",
    "    return functools.reduce(lambda a, b: a + b, ops) # Uses functools.reduce to sum these interaction terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "N5_C =  2*np.pi* np.array([[-0.36082693, -0.0702137 , -1.41518116],\n",
    "      [-0.0702137 , -0.60153649,  0.32312139],\n",
    "      [-1.41518116,  0.32312139, 50.80213093]]) # in MHz\n",
    "\t  \n",
    "N1_C = 2*np.pi*np.array([[  2.13814981,   3.19255832,  -2.48895215],\n",
    "      [  3.19255832,  15.45032887, -12.44778343],\n",
    "      [ -2.48895215, -12.44778343,  12.49532827]]) # in MHz\n",
    "\n",
    "N5_D =  2*np.pi*np.array([[-2.94412424e-01, -5.68059200e-02, -1.02860888e+00],\n",
    "      [-5.68059200e-02, -5.40578469e-01, -2.67686240e-02],\n",
    "      [-1.02860888e+00, -2.67686240e-02,  5.05815320e+01]]) # in MHz\n",
    "\t  \n",
    "N1_D = 2*np.pi* np.array([[ 0.98491908,  3.28010265, -0.53784491],\n",
    "      [ 3.28010265, 25.88547678, -1.6335986 ],\n",
    "      [-0.53784491, -1.6335986 ,  1.41368001]]) # in MHz\n",
    "\n",
    "\n",
    "ErC_Dee =  np.array([[ 26.47042689, -55.90357828,  50.1679204 ],\n",
    "                            [-55.90357828, -20.86385225,  76.13493805],\n",
    "                             [ 50.1679204,  76.13493805,  -5.60657464]]) # in Mrad/s\n",
    "\n",
    "\n",
    "\n",
    "ErD_Dee = np.array([[ 11.08087889, -34.6687169,   12.14623706],\n",
    "                            [-34.6687169,  -33.09039672,  22.36229081],\n",
    "                            [ 12.14623706,  22.36229081,  22.00951783]]) #  in Mrad/s\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "b0 = 1.4 * 2*math.pi # Zeeman field strength in radians per microsecond\n",
    "\n",
    "kr = 1. # Rate constant 1/us\n",
    "tmax = 12. / kr # Maximum time us \n",
    "tlist = np.linspace(0, tmax, math.ceil(1000*tmax)) # Time points for simulation\n",
    "theta = np.linspace(0, np.pi, 15)\n",
    "\n",
    "xyz = []\n",
    "\n",
    "for theta_ in theta:\n",
    "    x = np.sin(theta_)\n",
    "    y = 0\n",
    "    z = np.cos(theta_)\n",
    "    xyz.append([x, y, z])       \n",
    "\n",
    "oris = np.asarray(xyz)\n",
    "\n",
    "# Initialize arrays for latitude and longitude\n",
    "num_points = len(oris)\n",
    "lat = np.zeros(num_points)\n",
    "lon = np.zeros(num_points)\n",
    "\n",
    "# Convert Cartesian coordinates to latitude and longitude\n",
    "for i in range(num_points):\n",
    "    x, y, z = oris[i]\n",
    "    lat[i] = np.arcsin(z) * (180/np.pi)\n",
    "    lon[i] = np.arctan2(y, x) * (180/np.pi)\n",
    "\n",
    "dims = [2, 2, 3, 3, 3] # Dimensions of the system components (2 qubits, 1 spin-1 nucleus)\n",
    "dim = np.prod(dims) # Total dimension of the composite system\n",
    "\n",
    "yields  = []\n",
    "\n",
    "# Now H_C_list and H_D_list contain Hamiltonians for each orientation\n",
    "\n",
    "Ps = 1/4 * mkSpinOp(dims,[]) - mkH12(dims, 0, 1, np.identity(3)) # Singlet projection operator\n",
    "\n",
    "\n",
    "rho0_C = (Ps / Ps.tr()).full().flatten()# Initial density matrix, normalized projection operator for the singlet state.\n",
    "rho0_D = np.zeros_like(rho0_C)\n",
    "Ps = Ps.data\n",
    "\n",
    "# Combine the initial states into one vector\n",
    "initial_state = np.concatenate((rho0_C, rho0_D))\n",
    "\n",
    "sol = []\n",
    "\n",
    "for ori in oris:\n",
    "    B0 = b0 * ori  # Magnetic field vector along orientation\n",
    "    \n",
    "    # Compute Hamiltonians for each orientation\n",
    "    Hzee = mkH1(dims, 0, B0) + mkH1(dims, 1, B0) # Zeeman Hamiltonian for two spins\n",
    "    Hhfc_C = mkH12(dims, 0, 2, N5_C) + mkH12(dims, 1, 3, N1_C)\n",
    "    Hhfc_D = mkH12(dims, 0, 2, N5_D) + mkH12(dims, 1, 4, N1_D)\n",
    "    Hdee_C = mkH12(dims, 0, 1, ErC_Dee)\n",
    "    Hdee_D = mkH12(dims, 0, 1, ErD_Dee)\n",
    "    H0_C = Hzee + Hhfc_C + Hdee_C  # Total Hamiltonian for component C\n",
    "    H0_D = Hzee + Hhfc_D + Hdee_D  # Total Hamiltonian for component D\n",
    "    LeffC = -1j*qt.spre(H0_C) + 1j*qt.spost(H0_C.conj().trans())\n",
    "    LeffD = -1j*qt.spre(H0_D) + 1j*qt.spost(H0_D.conj().trans())\n",
    "    \n",
    "    mat = qt.Qobj(np.block([[LeffC.data - np.eye(LeffC.shape[0]), np.eye(LeffC.shape[0])], [np.eye(LeffD.shape[0]), LeffD.data - np.eye(LeffD.shape[0])]]))\n",
    "\n",
    "    sols = sci.sparse.linalg.spsolve(mat.data, initial_state)\n",
    "    sol = sol.append(sols)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}