Add files via upload

Author: CSQC-Duke

CR_CoherentMixing.ipynb

@@ -0,0 +1 @@
+{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[]},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["###**Model calculation: Idealized double-slit experiment**\n","\n","\n","\n","> This code can be used to reproduce the data plotted in Figure 1 of Schultz et al., Coherence in chemistry: foundations and frontiers. *Chem. Revs.* (submitted)\n","\n","To begin, press the play button beside each cell (proceed sequentially through the cells)."],"metadata":{"id":"H2WPl-rvZSYR"}},{"cell_type":"code","source":["#@title Import libraries, prep plot settings\n","\n","import matplotlib.pyplot as plt\n","from scipy import signal\n","import numpy as np\n","import math\n","from IPython.display import set_matplotlib_formats\n","set_matplotlib_formats('svg')\n","\n","#packages for plotting\n","import matplotlib\n","import matplotlib.pyplot as plt"],"metadata":{"id":"-X1cTwAH1cKg","cellView":"form"},"execution_count":null,"outputs":[]},{"cell_type":"code","metadata":{"id":"Nu_OWXfHeMFE","executionInfo":{"status":"ok","timestamp":1692576188156,"user_tz":240,"elapsed":197,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"source":["#@title Tunable parameters\n","\n","wavelength = 4.0\n","d = 20\n","L = 7"],"execution_count":2,"outputs":[]},{"cell_type":"code","metadata":{"cellView":"form","id":"nrAeY-12ZdST"},"source":["#@title Calculate and plot\n","\n","def Amplitude(L, y, d, wavelength): #calculates the AMPLITUDE of the sum of two sine waves y1 and y2 with amplitude a1 and a2, identical frequencies, and phases alpha1 and alpha2, respectively. For derivation, see Jenkins and White, 4th ed., pp 239\n","    R_1 = math.sqrt(L**2 + (y-d/2)**2)\n","    R_2 = math.sqrt(L**2 + (y+d/2)**2)\n","    a1 = 1/R_1 # INTENSITY goes as 1/R^2, not the amplitude\n","    a2 = 1/R_2\n","    alpha1 = 2*math.pi*R_1/wavelength\n","    alpha2 = 2*math.pi*R_2/wavelength\n","    sumamplitude = math.sqrt(a1**2 + a2**2 + 2*a1*a2*math.cos(alpha1 - alpha2))\n","    return sumamplitude\n","\n","y = np.arange(- 30, 30, 0.01)\n","\n","data = [Amplitude(L, y, d, wavelength)**2 for y in y]\n","data2 = [1/(L**2 + (y - d/2)**2) for y in y]\n","data3 = [1/(L**2 + (y + d/2)**2) for y in y]\n","data4 = [sum(x) for x in zip(data2, data3)]\n","\n","\n","fig = plt.figure(figsize=(6, 3))\n","ax = fig.add_subplot(1, 1, 1)\n","plt.plot(y,data, 'black', y, data4, 'b--')\n","plt.axis('off')\n","plt.show()\n","plt.savefig(\"Interference.svg\")"],"execution_count":null,"outputs":[]}]}

CR_DoubleSlit.ipynb

@@ -0,0 +1 @@
+{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[]},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["###**Model calculation: Interferometer**\n","\n","\n","\n","> This code can be used to reproduce the data plotted in Figure 2 of Schultz et al., Coherence in chemistry: foundations and frontiers. *Chem. Revs.* (submitted)\n","\n","To begin, press the play button beside each cell (proceed sequentially through the cells)."],"metadata":{"id":"H2WPl-rvZSYR"}},{"cell_type":"code","source":["#@title Import libraries, prep plot settings\n","\n","import numpy as np\n","import math\n","import matplotlib.pyplot as plt\n","from scipy import signal\n","from IPython.display import set_matplotlib_formats\n","set_matplotlib_formats('svg')\n","\n","def gensignal(N, trange, Autcrr):\n","    alpha_real = np.random.normal(0, 1/4, int(N))\n","    alpha_im = np.random.normal(0,1/4, int(N))\n","\n","    #initialize empty complex array to store random gaussian samples\n","    alpha = np.zeros(N+1)\n","    alpha = alpha.astype(complex)\n","\n","    for mu in range(1, int(N/2)+1):\n","        real = alpha_real[mu]\n","        im = 1j*alpha_im[mu]\n","        alpha[mu] = real + im\n","        alpha[-mu] = np.conjugate(alpha[mu])\n","\n","    # create alpha for zero freq component\n","    alpha[0] = 1\n","\n","    t= np.linspace(0, trange, N+1)\n","    #construct the spectral representation of the signal A using the random samples alpha_w\n","    fourtransAutcrr = np.fft.fft(Autcrr)\n","    A_w = [alpha[mu]*fourtransAutcrr[mu] for mu in range(-int(N/2), int(N/2)+1)]\n","\n","\n","\n","\n","    A_w = np.fft.ifftshift(A_w)\n","\n","    A_t = np.fft.ifft(A_w)\n","    return np.real(A_t) #discard small imaginary part\n","\n","#calculates the autocorrelation function signal[i]*signal[i+j] averaged over Neave different i's\n","def autocorrelation(signal, j, Neave):\n","    sum = 0\n","    signalj = np.roll(signal, -j)       #creates a copy of the timeseries and shifts it over by j\n","    signal = signal[:Neave]                #remove data from lists that are not included in the ergotic average (Neave is number of timestamps in ergotic average)\n","    signalj= signalj[:Neave]\n","    return np.dot(signal, signalj)/Neave    #returns the ergotic average of signal[i]*signal[i+j] over Neave different i's"],"metadata":{"id":"IksVGVrZzYus","cellView":"form"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Tunable parameters\n","\n","N=2**20\n","dt = 0.01\n","trange = (N+1)*dt\n","t= np.linspace(0, trange, N+1)\n","\n","epsilon = 20\n","tau =   20\n","period = 1\n","omega = 2*math.pi/period"],"metadata":{"id":"-X1cTwAH1cKg"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Calculate and plot\n","\n","\n","###############\n","## Calculate main pulse\n","###############\n","Autcrr = [np.cos(omega*t)*epsilon/tau*np.exp(-t/tau) for t in t]\n","sig = gensignal(N, trange, Autcrr)\n","\n","\n","## Set up delay vectors\n","###################\n","Ndelaysamples = 150 #number of data points in delay plot\n","tdelayrange = 8*period\n","dtdelay = tdelayrange/Ndelaysamples\n","dN = int(dtdelay/dt)\n","t_delay = np.linspace(0, tdelayrange, Ndelaysamples)\n","\n","\n","integrationtime = 2000 * tau\n","Ninteg = int(integrationtime/dt)\n","\n","\n","#Calculate the integrated power on the detector sig[i]sig[i+j] where j is the delay and i sums over the integration time\n","##########\n","detectorsig = np.zeros(Ndelaysamples)\n","\n","for delay in range(0, Ndelaysamples):\n","    sigdelay = np.roll(sig, -delay*dN)\n","    sigdelaydetect = sigdelay[:Ninteg]\n","    sigdetect = sig[:Ninteg]\n","    intensitysum = [(sigdelaydetect[j]+sigdetect[j])**2 for j in range(len(sigdetect))] #instantaneous intensity at time j recorded over the integration time\n","    detectorsig[delay] = np.mean(intensitysum) #average intensity\n","\n","\n","intensitytest = [sig[j]**2 for j in range(len(sig))]\n","val = 2*np.mean(intensitytest)\n","incoherentsum = [val for j in range(Ndelaysamples)]\n","\n","\n","############\n","## Plotting\n","############\n","\n","fig = plt.figure(figsize = (3,6))\n","\n","ax1 = fig.add_subplot(3,1,1)\n","ax1.plot(t, sig, linewidth = 1, color = 'b', label = 'Pulse')\n","ax1.plot(t, np.roll(sig, int(2.4*period/dt)), linewidth = 1, color = 'g', label = 'Delay') #comparison between reference and sample delay\n","plt.xlim(0, 7*period)\n","ax1.axes.xaxis.set_ticks([])\n","ax1.axes.yaxis.set_ticks([])\n","ax1.set_xlabel('Time')\n","ax1.set_ylabel('Signal')\n","ax1.set_title('Pulse Profiles')\n","ax1.legend(prop={'size': 7})\n","\n","ax2 = fig.add_subplot(3,1,2)\n","ax2.plot(t_delay, detectorsig, linewidth = 1, color = 'black', label = 'Interference')\n","ax2.plot(t_delay, incoherentsum, linewidth = 1.5, color = 'blue', linestyle = ':', label = 'Incoherent sum')\n","ax2.set_xlabel('Delay time')\n","ax2.set_ylabel('Average Intensity')\n","ax2.axes.xaxis.set_ticks([])\n","ax2.axes.yaxis.set_ticks([])\n","ax2.legend(prop={'size': 8})\n","ax2.set_title('Detector response')\n","fig.tight_layout()\n","\n","\n","fig2 =  plt.figure(figsize = (1,1))\n","ax3 = fig2.add_subplot(1,1,1)\n","x = np.linspace(0,1, int((N+1)/2)+1)\n","y = np.abs(np.fft.rfft(sig))\n","ax3.plot(x, y, linewidth = 1, color = 'black')\n","plt.xlim(0, 0.05)\n","ax3.axes.xaxis.set_ticks([])\n","ax3.axes.yaxis.set_ticks([])\n","ax3.set_xlabel('Frequency',fontdict={'fontsize': 8})\n","ax3.set_ylabel('Power',fontdict={'fontsize': 8})"],"metadata":{"id":"f932zxF_Wx0b","cellView":"form"},"execution_count":null,"outputs":[]}]}

CR_Interference.ipynb

@@ -0,0 +1 @@
+{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[{"file_id":"1kTrDRr1Logxh8_dirw7fxJkfeOXs_WWW","timestamp":1644613029419}]},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["###**Model calculation: Molecular interferometer**\n","  \n","> This code can be used to reproduce the data plotted in Figure 21 of Schultz et al., Coherence in chemistry: foundations and frontiers. *Chem. Revs.* (submitted)\n","\n","To begin, press the play button beside each cell (proceed sequentially through the cells)."],"metadata":{"id":"H2WPl-rvZSYR"}},{"cell_type":"code","source":["#@title Import libraries, prep plot settings\n","\n","import matplotlib.pyplot as plt\n","import numpy as np\n","from math import *\n","import cmath\n","from scipy import linalg\n","from IPython.display import set_matplotlib_formats\n","set_matplotlib_formats('svg')\n","\n","hbar = 1\n","\n","i_ = complex(0, 1)\n","\n","SMALL_SIZE = 18\n","MEDIUM_SIZE = 20\n","BIGGER_SIZE = 22\n","\n","plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n","plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n","plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n","plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n","plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n","plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n","plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title"],"metadata":{"id":"-X1cTwAH1cKg","cellView":"form"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Tunable parameters\n","\n","V = 0.5\n","V_phase_L = 0\n","V_phase_R = np.arange(0,4*pi,0.1,dtype=float)\n","site_energies = [0, 2, 2, -2]\n","points = 1000"],"metadata":{"id":"FDhVIEa-NHzX","executionInfo":{"status":"ok","timestamp":1692631038460,"user_tz":240,"elapsed":9,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"execution_count":95,"outputs":[]},{"cell_type":"code","source":["#@title Set up arrays and functions\n","Delta_phi = V_phase_R-V_phase_L\n","\n","DB1coupling = V*np.exp(-i_*V_phase_L)\n","B1Acoupling = V*np.exp(-i_*V_phase_L)\n","B2Acoupling = V*np.exp(-i_*V_phase_L)\n","\n","rho_init = np.array([[1, 0, 0, 0], [0, 0, 0, 0], [ 0, 0, 0, 0], [0, 0, 0, 0]], dtype='complex') # Initial density matrix\n","sys_size = rho_init.shape[0]\n","\n","def dag(A):\n","  return np.conjugate(A)\n","\n","def observable(P,B):\n","  return np.real(np.trace(B @ P))    # expectation value is the trace of the projection (P) onto B (rho eig)\n","\n","proj = np.zeros([sys_size,sys_size])  # Projection onto site 4 (or the population on the acceptor) is the observable\n","proj[3,3] = 1"],"metadata":{"id":"_mn6Z6jdY_zE","cellView":"form","executionInfo":{"status":"ok","timestamp":1692631040717,"user_tz":240,"elapsed":188,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"execution_count":96,"outputs":[]},{"cell_type":"code","source":["#@title Calculate unitary time-evolution of density matrix (in the site basis)\n","\n","maxpop = []\n","for n in range(len(V_phase_R)):\n","\n","  DB2coupling = V*np.exp(-i_*V_phase_R[n])\n","\n","  H = np.array([[site_energies[0], DB1coupling, DB2coupling, 0], [dag(DB1coupling), site_energies[1], 0, B1Acoupling], [dag(DB2coupling), 0, site_energies[2], B2Acoupling], [0, dag(B1Acoupling), dag(B2Acoupling), site_energies[3]]], dtype='complex')   # Hamiltonian\n","\n","  maxEdiff = np.max(np.abs(np.diff(H)))\n","  tmax = 25/maxEdiff\n","  t = np.linspace(0,tmax,points)    # Time vector\n","  dt = t[1]-t[0]                    # Time step\n","\n","  U = linalg.expm(-i_*H*dt)         # Propagators\n","  Udag = linalg.expm(i_*H*dt)\n","\n","  rho_store = []\n","  rho_store.append(rho_init)\n","  population = []\n","  population.append(observable(proj,rho_store[0]))\n","  rho_time = rho_init.copy()\n","\n","  for i in range(len(t)-1):\n","    k = i+1\n","    rho_time = U @ rho_time @ Udag                      # Propagate the density matrix by dt\n","    rho_store.append(rho_time)                          # Store propagated density matrix\n","    population.append(observable(proj,rho_store[k]))    # Calculate the observable\n","\n","  maxpop.append(np.max(population))"],"metadata":{"id":"LM5_ZsZnat77","cellView":"form","executionInfo":{"status":"ok","timestamp":1692631049501,"user_tz":240,"elapsed":5777,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"execution_count":97,"outputs":[]},{"cell_type":"code","source":["#@title Plot\n","\n","fig = plt.figure(figsize=(15,6))\n","ax1 = fig.add_subplot(1, 2, 1)\n","ax1.plot(t,population, 'r-', linewidth=1)\n","plt.title(\"Example population dynamics\")\n","plt.xlabel(\"time (arb. units)\")\n","plt.ylabel(\"Acceptor population\")\n","\n","ax2 = fig.add_subplot(1, 2, 2)\n","ax2.plot(Delta_phi,maxpop, 'r-', linewidth=1)\n","ax2.xaxis.set_ticks([0,pi,2*pi,3*pi, 4*pi])\n","plt.title(\"Coupling pathway interference\")\n","plt.xlabel(\"Phase difference ($\\Delta \\phi_{LR}$)\")\n","plt.ylabel(\"Acceptor population\")\n","\n","fig.tight_layout()"],"metadata":{"id":"nonR36BmM4aE","cellView":"form"},"execution_count":null,"outputs":[]}]}

CR_MolecularInterferometer.ipynb

@@ -0,0 +1 @@
+{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[{"file_id":"1kTrDRr1Logxh8_dirw7fxJkfeOXs_WWW","timestamp":1644613029419}]},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","source":["###**Model calculation: Rabi oscillations and quantum beats**\n","\n","\n","\n","> This code can be used to reproduce the data plotted in Figures 10 and 11 of Schultz et al., Coherence in chemistry: foundations and frontiers. *Chem. Revs.* (submitted)\n","\n","To begin, press the play button beside each cell (proceed sequentially through the cells). *Note: run either the \"site 1\" or \"eigenstate 1\" initial state cells and compare the calculated time evolution.*"],"metadata":{"id":"H2WPl-rvZSYR"}},{"cell_type":"code","source":["#@title Import libraries, prep plot settings\n","\n","import matplotlib.pyplot as plt\n","import numpy as np\n","from math import *\n","import cmath\n","from scipy import linalg\n","from IPython.display import set_matplotlib_formats\n","set_matplotlib_formats('svg')\n","\n","hbar = 1\n","\n","i_ = complex(0, 1)\n","\n","SMALL_SIZE = 18\n","MEDIUM_SIZE = 20\n","BIGGER_SIZE = 22\n","\n","plt.rc('font', size=SMALL_SIZE)          # controls default text sizes\n","plt.rc('axes', titlesize=SMALL_SIZE)     # fontsize of the axes title\n","plt.rc('axes', labelsize=MEDIUM_SIZE)    # fontsize of the x and y labels\n","plt.rc('xtick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n","plt.rc('ytick', labelsize=SMALL_SIZE)    # fontsize of the tick labels\n","plt.rc('legend', fontsize=SMALL_SIZE)    # legend fontsize\n","plt.rc('figure', titlesize=BIGGER_SIZE)  # fontsize of the figure title"],"metadata":{"id":"Cv5auRqrYQme","cellView":"form"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Tunable parameters\n","\n","n = 1000                 # number of points\n","tmax = 8                 # max time value\n","E1 = 10\n","E2 = 10\n","V = 0.05*E1"],"metadata":{"id":"bSUnrx6SE5rl","executionInfo":{"status":"ok","timestamp":1692631157079,"user_tz":240,"elapsed":167,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"execution_count":2,"outputs":[]},{"cell_type":"code","source":["#@title Set up arrays\n","\n","dt = tmax/n                # timestep\n","t = np.arange(0,tmax,dt)   # time vector\n","\n","H = np.array([[E1, V], [V, E2]], dtype='complex')\n","\n","eigval, eigvec = linalg.eigh(H)"],"metadata":{"id":"_mn6Z6jdY_zE","cellView":"form","executionInfo":{"status":"ok","timestamp":1692631158422,"user_tz":240,"elapsed":149,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"execution_count":3,"outputs":[]},{"cell_type":"code","source":["#@title Initial state = site 1, projection to sites\n","\n","rho = np.array([[1, 0], [0, 0]], dtype='complex')\n","rho_eig = eigvec @ rho @ np.conjugate(eigvec.transpose())\n","\n","Proj1 = np.array([[1, 0], [0, 0]], dtype='complex')\n","Proj2 = np.array([[0, 0], [0, 1]], dtype='complex')\n","Proj1_eig = eigvec @ Proj1 @ np.conjugate(eigvec.transpose())\n","Proj2_eig = eigvec @ Proj2 @ np.conjugate(eigvec.transpose())"],"metadata":{"id":"P4v02bBtstsy","cellView":"form","executionInfo":{"status":"ok","timestamp":1692631161331,"user_tz":240,"elapsed":197,"user":{"displayName":"Quantum Coherence","userId":"04118399231983577771"}}},"execution_count":4,"outputs":[]},{"cell_type":"code","source":["#@title Initial state = eigenstate 1, projection to sites\n","\n","rho_eig = np.array([[1, 0], [0, 0]], dtype='complex')\n","\n","Proj1 = np.array([[1, 0], [0, 0]], dtype='complex')\n","Proj2 = np.array([[0, 0], [0, 1]], dtype='complex')\n","Proj1_eig = eigvec @ Proj1 @ np.conjugate(eigvec.transpose())\n","Proj2_eig = eigvec @ Proj2 @ np.conjugate(eigvec.transpose())"],"metadata":{"cellView":"form","id":"dp8TC20_szLz"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Calculate UNITARY time-evolution of density matrix\n","\n","def U(tn):\n","  return linalg.expm(-i_*np.diag(eigval)*tn)\n","\n","def Udag(tn):\n","  return np.transpose(np.conjugate(linalg.expm(-i_*np.diag(eigval)*tn)))\n","\n","P1 = np.zeros((n,))\n","P2 = np.zeros((n,))\n","rho_store = []\n","test = np.zeros((n,))\n","\n","for i in range(n):\n","   rho_time = U(t[i]) @ rho_eig.copy() @ Udag(t[i])     # Propagate the density matrix\n","   P1[i] = np.real(np.trace(rho_time @ Proj1_eig))  # Calculate the observable at time t\n","   P2[i] = np.real(np.trace(rho_time @ Proj2_eig))  # Calculate the observable at time t\n","   rho_store.append(rho_time)\n","\n","fig = plt.figure(figsize=(9,6))\n","ax = fig.add_subplot(1, 1, 1)\n","ax.plot(t,np.real(P1), 'r-', linewidth=1)\n","ax.plot(t,np.real(P2), 'r--', linewidth=3)\n","ax.xaxis.set_ticks([0])\n","ax.yaxis.set_ticks([0,1])\n","plt.xlabel(\"Time (arb. units)\")\n","plt.ylabel(\"Population\")\n","plt.legend(['$\\phi_1$', '$\\phi_2$'])\n","fig.tight_layout()\n","plt.savefig(\"output.svg\")"],"metadata":{"id":"LM5_ZsZnat77","cellView":"form"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["#@title Calculate beats for driven 3LS\n","\n","H3LS = np.array([[0, 0, 0], [0, E1, V], [0, V, E2]], dtype='complex')\n","eigval, eigvec = linalg.eigh(H3LS)\n","\n","total_signal = 2*np.cos((eigval[2]+eigval[1])*t)*np.cos(abs(eigval[2]-eigval[1])*t)\n","average_signal = 2*np.cos((eigval[2]+eigval[1])*t)\n","beat_signal = 2*np.cos(abs(eigval[2]-eigval[1])*t)\n","\n","\n","fig = plt.figure(figsize=(8,5))\n","ax = fig.add_subplot(1, 1, 1)\n","ax.plot(t,average_signal, 'b--', label =\"average\", linewidth=0.5)\n","ax.plot(t,beat_signal, 'r', label =\"beat\", linewidth=4)\n","ax.plot(t,total_signal, 'k', label =\"detected\", linewidth=2)\n","ax.xaxis.set_ticks([0])\n","ax.yaxis.set_ticks([])\n","plt.xlabel(\"Time (arb. units)\")\n","plt.ylabel(\"I(t) (arb. units)\")\n","plt.legend(['total', 'mean', 'beat'],loc =\"right\")\n","fig.tight_layout()\n","plt.savefig(\"output.svg\")"],"metadata":{"id":"jlftZRbX2nVW","cellView":"form"},"execution_count":null,"outputs":[]}]}