Simulate three-wave mixing effect for chirality detection using data from Nature 497 (2013) 475
Here we recycle codes from the previous exercises “Compute rotational energy levels of rigid water molecule” and “Compute Stark effect for rotational levels of water molecule” and in the first step adapt them to a chiral molecule 1,2-propanediol
[16]:
# Take molecular data for 1,2-propanediol from Nature 497 (2013) 475
# rotational constants in MHz
A = 8572.05
B = 3640.10
C = 2790.96
# is the molecule near-prolate or near-oblate?
# dipole moment in the principal axes system in units of Debye
dipole_R = [1.2, 1.9, 0.36] # R-enantiomer
dipole_S = [1.2, 1.9, -0.36] # S-enantiomer
The following part is copied from the exercise “Compute rotational energy levels of rigid water molecule”
[17]:
import numpy as np
# Set up functions for operators J+|psi>, J-|psi>, Jz|psi>, J^2|psi> and overlap
delta = lambda a, b: 1.0 if a==b else 0.0
Jplus = lambda J, k, c=1: (J, k-1, np.sqrt(J * (J + 1) - k * (k - 1)) * c if abs(k-1) <= J else 0)
Jminus = lambda J, k, c=1: (J, k+1, np.sqrt(J * (J + 1) - k * (k + 1)) * c if abs(k+1) <= J else 0)
Jz = lambda J, k, c=1: (J, k, k * c)
Jsq = lambda J, k, c=1: (J, k, J * (J + 1) * c)
# function for overlap integral <psi'|operator psi>
overlap = lambda J1, k1, c1, J2, k2, c2: c1 * c2 * delta(J1, J2) * delta(k1, k2)
# Compute rotaitonal energies and wave functions
Jmax = 10 # maximal value of J quantum number spanned by basis
enr = {} # stores energies for different J quanta
vec = {} # stores eigenvectors for different J quanta
k_quanta = {} # stores list of k quantun numbers for different J quanta
for J in range(Jmax+1):
# generate basis set indices (k-quanta)
k_list = [k for k in range(-J, J+1)]
# matrix elements of J+^2, J_^2, Jz^2, and J^2
Jplus_mat = np.array([[overlap(*(J, k1, 1), *Jplus(*Jplus(J, k2))) for k1 in k_list] for k2 in k_list], dtype=np.float64)
Jminus_mat = np.array([[overlap(*(J, k1, 1), *Jminus(*Jminus(J, k2))) for k1 in k_list] for k2 in k_list], dtype=np.float64)
Jz_mat = np.array([[overlap(*(J, k1, 1), *Jz(*Jz(J, k2))) for k1 in k_list] for k2 in k_list], dtype=np.float64)
Jsq_mat = np.array([[overlap(*(J, k1, 1), *Jsq(J, k2)) for k1 in k_list] for k2 in k_list], dtype=np.float64)
# Hamiltonian for near-prolate top
hmat_a = (Jplus_mat + Jminus_mat) * (B - C)/4.0 + Jz_mat * (2*A - B - C)/2.0 + (B + C)/2.0 * Jsq_mat
enr_a, vec_a = np.linalg.eigh(hmat_a)
# Hamiltonian for near-oblate top
hmat_c = (Jplus_mat + Jminus_mat) * (A - B)/4.0 + Jz_mat * (2*C - A - B)/2.0 + (A + B)/2.0 * Jsq_mat
enr_c, vec_c = np.linalg.eigh(hmat_c)
# print energies and assignments by k_a and k_c quantum numbers
for istate in range(len(k_list)):
ind_a = np.argmax(vec_a[:,istate]**2)
ind_c = np.argmax(vec_c[:,istate]**2)
print(J, istate, " %12.6f"%enr_a[istate], "ka, kc = ", abs(k_list[ind_a]), abs(k_list[ind_c]))
enr[J] = enr_a
vec[J] = vec_a
k_quanta[J] = k_list
0 0 0.000000 ka, kc = 0 0
1 0 6431.060000 ka, kc = 0 1
1 1 11363.010000 ka, kc = 1 1
1 2 12212.150000 ka, kc = 1 0
2 0 19192.694106 ka, kc = 0 2
2 1 23375.990000 ka, kc = 1 2
2 2 25923.410000 ka, kc = 1 1
2 3 40719.260000 ka, kc = 2 1
2 4 40819.745894 ka, kc = 2 0
3 0 38092.937260 ka, kc = 0 3
3 1 41335.980211 ka, kc = 1 3
3 2 46423.325008 ka, kc = 1 0
3 3 60012.440000 ka, kc = 2 2
3 4 60505.862740 ka, kc = 2 1
3 5 86854.519789 ka, kc = 3 1
3 6 86862.014992 ka, kc = 3 0
4 0 62886.983186 ka, kc = 0 4
4 1 65181.519337 ka, kc = 1 4
4 2 73620.870782 ka, kc = 1 3
4 3 85658.253452 ka, kc = 2 3
4 4 87081.407179 ka, kc = 2 0
4 5 112759.180663 ka, kc = 3 2
4 6 112811.229218 ka, kc = 3 1
4 7 150093.346548 ka, kc = 4 1
4 8 150093.809635 ka, kc = 4 0
5 0 93352.854217 ka, kc = 0 5
5 1 94845.390415 ka, kc = 1 5
5 2 107377.573123 ka, kc = 1 4
5 3 117590.522318 ka, kc = 2 4
5 4 120699.422272 ka, kc = 2 1
5 5 145188.269017 ka, kc = 3 3
5 6 145393.160453 ka, kc = 3 0
5 7 182471.677682 ka, kc = 4 2
5 8 182475.823511 ka, kc = 4 3
5 9 230473.690568 ka, kc = 5 1
5 10 230473.716425 ka, kc = 5 0
6 0 129357.316924 ka, kc = 0 6
6 1 130263.373854 ka, kc = 1 6
6 2 147497.709775 ka, kc = 1 5
6 3 155729.089305 ka, kc = 2 5
6 4 161403.491863 ka, kc = 2 0
6 5 184142.835841 ka, kc = 3 4
6 6 184740.156676 ka, kc = 3 5
6 7 221394.598687 ka, kc = 4 3
6 8 221415.137849 ka, kc = 4 4
6 9 269312.800305 ka, kc = 5 2
6 10 269313.083549 ka, kc = 5 3
6 11 327998.212009 ka, kc = 6 1
6 12 327998.213364 ka, kc = 6 0
7 0 170858.758031 ka, kc = 0 7
7 1 171380.755042 ka, kc = 1 7
7 2 193730.607446 ka, kc = 1 6
7 3 199983.211126 ka, kc = 2 6
7 4 209119.915561 ka, kc = 2 1
7 5 229601.930471 ka, kc = 3 5
7 6 231026.309741 ka, kc = 3 0
7 7 266892.143447 ka, kc = 4 4
7 8 266966.343428 ka, kc = 4 5
7 9 314688.791198 ka, kc = 5 3
7 10 314690.479457 ka, kc = 5 4
7 11 373298.805427 ka, kc = 6 2
7 12 373298.822980 ka, kc = 6 3
7 13 442666.963289 ka, kc = 7 1
7 14 442666.963357 ka, kc = 7 0
8 0 217866.001181 ka, kc = 0 8
8 1 218155.391660 ka, kc = 1 8
8 2 245801.514781 ka, kc = 1 7
8 3 250256.644928 ka, kc = 2 7
8 4 263691.090165 ka, kc = 2 2
8 5 281519.940136 ka, kc = 3 6
8 6 284435.612608 ka, kc = 3 1
8 7 318989.645820 ka, kc = 4 3
8 8 319207.237179 ka, kc = 4 6
8 9 366630.401464 ka, kc = 5 6
8 10 366637.644854 ka, kc = 5 5
8 11 425128.854590 ka, kc = 6 5
8 12 425128.976810 ka, kc = 6 4
8 13 494430.066741 ka, kc = 7 2
8 14 494430.067757 ka, kc = 7 3
8 15 574479.894663 ka, kc = 8 1
8 16 574479.894666 ka, kc = 8 0
9 0 270401.652883 ka, kc = 0 9
9 1 270557.572427 ka, kc = 1 9
9 2 303472.977068 ka, kc = 1 8
9 3 306454.001227 ka, kc = 2 8
9 4 324902.453597 ka, kc = 0 5
9 5 339827.712347 ka, kc = 3 7
9 6 345098.608081 ka, kc = 3 0
9 7 377701.332128 ka, kc = 4 4
9 8 378248.320325 ka, kc = 4 7
9 9 425168.977949 ka, kc = 5 7
9 10 425193.969478 ka, kc = 5 6
9 11 483511.717027 ka, kc = 6 6
9 12 483512.323521 ka, kc = 6 5
9 13 552717.410875 ka, kc = 7 5
9 14 552717.418971 ka, kc = 7 4
9 15 632706.149618 ka, kc = 8 2
9 16 632706.149674 ka, kc = 8 3
9 17 723436.976402 ka, kc = 9 1
9 18 723436.976403 ka, kc = 9 0
10 0 328485.690589 ka, kc = 0 10
10 1 328567.884822 ka, kc = 1 10
10 2 366600.815746 ka, kc = 1 9
10 3 368487.264844 ka, kc = 2 9
10 4 392491.870458 ka, kc = 0 6
10 5 404436.603892 ka, kc = 3 8
10 6 413032.834697 ka, kc = 3 1
10 7 443023.793888 ka, kc = 4 5
10 8 444239.396405 ka, kc = 4 8
10 9 490336.491392 ka, kc = 5 8
10 10 490409.984066 ka, kc = 5 7
10 11 548473.734838 ka, kc = 6 7
10 12 548476.135614 ka, kc = 6 0
10 13 617548.182885 ka, kc = 7 6
10 14 617548.228479 ka, kc = 7 5
10 15 697452.914711 ka, kc = 8 5
10 16 697452.915214 ka, kc = 8 4
10 17 788126.787009 ka, kc = 9 4
10 18 788126.787012 ka, kc = 9 3
10 19 889538.191719 ka, kc = 10 1
10 20 889538.191719 ka, kc = 10 0
At this point we have computed rotational states of 1,2-propanediol molecule for \(J=0..10\). In fact for simulations of the three-wave mixing effect we need only three states. Those used in the paper are (\(J\), \(k_a\), \(k_c\)) = (0, 0, 0), (1, 1, 0), and (1, 1, 1), which according to the printout above have the absolute indexes 0 (for \(J=0\)), 2 and 1 (for \(J=1\)). We are going to use only these states:
[18]:
# keep energies and eigenvectors of rotational states with indexes 0 (J=0), 2, and 1 (J=1)
enr3 = {0 : enr[0], 1 : enr[1][[2, 1]]} # keys=J, values=lists of k quanta
vec3 = {0 : vec[0], 1 : vec[1][:, [2, 1]]}
The following parts are copied from the exercise “Compute Stark effect for rotational levels of water molecule”
[19]:
# Set up functions for computing symmetric-top matrix elements of direction cosine matrix
import sympy
from sympy.physics.wigner import wigner_3j
# matrix of transformation from Cartesian to spherical-tensor form for rank-1 operator
# order of spherical components in rows: (1,1), (1,0), (1,-1)
# order of Cartesian components in columns: x, y, z
umat = np.array( [[-1.0/np.sqrt(2.0), -1j/np.sqrt(2.0), 0.0], \
[0.0, 0.0, 1.0], \
[1.0/np.sqrt(2.0), -1j/np.sqrt(2.0), 0.0]], dtype=np.complex128)
# split symmetric-top matrix element of direction consine matrix into a product of two parts M and K
# this is computationally much more efficient since we can now treat the m- and k-dependent parts (laboratory
# and molecular frames' quanta) of basis separately, while the complete basis is produced by a Kronecker
# product of all m and all k quanta, i.e., may become very large
# K-matrix contracted with permanent dipole moment, i.e. sum_alpha(K_alpha * dipole_alpha)
# R-enantiomer
def kmat_R(omega, J1, k1, J2, k2):
threej = np.array([wigner_3j(J2, omega, J1, k2, sigma, -k1) \
for sigma in reversed(range(-omega, omega + 1))], dtype=np.float64)
return (-1.0)**k1 * np.dot(np.dot(threej, umat), dipole_R)
# S-enantiomer
def kmat_S(omega, J1, k1, J2, k2):
threej = np.array([wigner_3j(J2, omega, J1, k2, sigma, -k1) \
for sigma in reversed(range(-omega, omega + 1))], dtype=np.float64)
return (-1.0)**k1 * np.dot(np.dot(threej, umat), dipole_S)
# M-matrix (does not depend on molecular parameters)
def mmat(omega, J1, m1, J2, m2):
threej = np.array([wigner_3j(J2, omega, J1, m2, sigma, -m1) \
for sigma in reversed(range(-omega, omega + 1))], dtype=np.float64)
return np.sqrt(2.0 * J1 + 1) * np.sqrt(2.0 * J2 + 1) * (-1.0)**m1 * np.dot(np.linalg.pinv(umat), threej)
[20]:
# Set up basis quanta and precompute matrices composing the total Hamiltonian
# the benefits of splitting the matrix element of direction cosine matrix becomes apparent here!
# list of J quanta spanned by basis
J_list = [J for J in enr3.keys()]
# list of m-quanta spanned by basis for each J
m_list = {J : [m for m in range(-J, J + 1)] for J in J_list}
# precompute matrix representation of Trot, which is diagonal in its own eigenbasis
h0 = {}
for J in J_list:
h0[J] = [e for m in m_list[J] for e in enr3[J]]
# precompute M and K matrices
mme = {}
kme_R = {}
kme_S = {}
for J1 in J_list:
for J2 in J_list:
# K-matrix in symmetric-top basis
me_R = np.array([[ kmat_R(1, J1, k1, J2, k2) for k2 in k_quanta[J2] ] for k1 in k_quanta[J1] ]) # shape = (k', k)
# transform K matrix to a basis of Trot eigenfunctions, i.e., rotational wavefunctions of water
kme_R[(J1, J2)] = np.dot(np.conjugate(vec3[J1].T), np.dot(me_R, vec3[J2])) # shape = (l', l)
# K-matrix in symmetric-top basis
me_S = np.array([[ kmat_S(1, J1, k1, J2, k2) for k2 in k_quanta[J2] ] for k1 in k_quanta[J1] ]) # shape = (k', k)
# transform K matrix to a basis of Trot eigenfunctions, i.e., rotational wavefunctions of water
kme_S[(J1, J2)] = np.dot(np.conjugate(vec3[J1].T), np.dot(me_S, vec3[J2])) # shape = (l', l)
# M-matrix in symmetric-top basis
mme[(J1, J2)] = np.array([[ mmat(1, J1, m1, J2, m2) for m2 in m_list[J2] ] for m1 in m_list[J1] ]) # shape = (m', m, A)
# matrix elements of dipole moment: dipole_A = M_A * K (where A = X, Y, Z)
# R-enantiomer
dipme_R = [ np.block([[ np.kron( mme[(J1, J2)][:,:,A], kme_R[(J1, J2)] ) \
for J2 in J_list] for J1 in J_list]) \
for A in range(3) ]
# S-enantiomer
dipme_S = [ np.block([[ np.kron( mme[(J1, J2)][:,:,A], kme_S[(J1, J2)] ) \
for J2 in J_list] for J1 in J_list]) \
for A in range(3) ]
Now we are ready set up the external microwave fields and the total Hamiltonian, and solve the TDSE
[30]:
from scipy import linalg
import sys
conv_to_MHz = 0.5034485879043021 # converts dipole[Debye] * field[V/cm] into MHz (CHECK!!!)
# note on unit conversion: Debye = 1/c * 10^{-21} Coulomb * meter and Coulomb * Volt = Joule
# microwave fields in V/cm
# along X - static field
field_x = lambda time: 65.0 if time < 800.0 else 0.0
# along Z - resonant Pi/2 pulse for 0->1 transition
omega = enr3[1][0] - enr3[0][0] # in MHz (= 1/s * 1e6 = 1/ps * 1e-12 * 1e6)
field_z = lambda time: 500*np.cos(2.0*np.pi*omega*1e-6*time) if time < 900 else 0
# total dimension of basis
dimen = sum(len(h0[J]) for J in J_list)
# time step and time grid in picoseconds
dt = 1
time_grid = np.linspace(0, 2000, 2000/dt)
# initial vector in the ground rotational state
wavepacket_R = np.zeros((dimen, len(time_grid)+1), dtype=np.complex128)
wavepacket_R[0,0] = 1.0
wavepacket_S = np.zeros((dimen, len(time_grid)+1), dtype=np.complex128)
wavepacket_S[0,0] = 1.0
# expectation values of dipole moment
dipmom_R = np.zeros((3, len(time_grid)+1), dtype=np.complex128)
dipmom_S = np.zeros((3, len(time_grid)+1), dtype=np.complex128)
# start propagation
for it, t in enumerate(time_grid):
field = [field_x(t), 0, field_z(t)] # field in V/cm
# H' = -dipole * field = -sum_A(M_A * F_A) * K, convert to units MHz
hmat_R = -np.block([[ np.kron( np.dot(mme[(J1, J2)], field), kme_R[(J1, J2)] ) * conv_to_MHz \
for J2 in J_list] for J1 in J_list])
hmat_S = -np.block([[ np.kron( np.dot(mme[(J1, J2)], field), kme_S[(J1, J2)] ) * conv_to_MHz \
for J2 in J_list] for J1 in J_list])
# H = H0 + H'
hmat_R += np.diag(np.concatenate([h0[J] for J in J_list]))
hmat_S += np.diag(np.concatenate([h0[J] for J in J_list]))
# time-evolution operator
expm_R = linalg.expm(-1j*dt*hmat_R*2.0*np.pi*1e-6)
expm_S = linalg.expm(-1j*dt*hmat_S*2.0*np.pi*1e-6)
# update wavepacket
wavepacket_R[:,it+1] = np.dot(expm_R, wavepacket_R[:,it])
wavepacket_S[:,it+1] = np.dot(expm_S, wavepacket_S[:,it])
# expectation values of dipole
dipmom_R[:,it+1] = np.dot(np.dot(dipme_R, wavepacket_R[:,it+1]), np.conjugate(wavepacket_R[:,it+1]))
dipmom_S[:,it+1] = np.dot(np.dot(dipme_S, wavepacket_S[:,it+1]), np.conjugate(wavepacket_S[:,it+1]))
[46]:
import matplotlib.pyplot as plt
ntimes = len(time_grid)
#plt.plot(time_grid, wavepacket_R[0,:ntimes].imag)
#plt.plot(time_grid, wavepacket_R[1,:ntimes].imag)
#plt.plot(time_grid, wavepacket_R[2,:ntimes].imag)
#plt.plot(time_grid, wavepacket_R[3,:ntimes].imag)
plt.plot(time_grid, dipmom_R[1,:ntimes])
plt.plot(time_grid, dipmom_S[1,:ntimes])
[46]:
[<matplotlib.lines.Line2D at 0x7f9b3c58f470>]
[ ]: