Compute rotational energy levels of rigid water molecule
Use the following geometrical parameters: the O–H distances are equal to 0.96 Angstrom and the H–O–H valence bond angle is 104\(^\circ\)
[13]:
# Compute Cartesain coordinates of atoms, moment of itertia, and rotational constants
import numpy as np
r1 = 0.96
r2 = 0.96
theta = 104.0 * np.pi / 180.0 # angle in radians
masses = [15.99491462, 1.00782503, 1.00782503] # masses of O, H1, and H2 atoms
# Cartesian coordinates of atoms (iatom, ix) where ix=0,1,2 (=x,y,z)
xyz = np.array( [[0, 0, 0], \
[r1 * np.sin(theta / 2.0), 0, r1 * np.cos(theta / 2.0)], \
[-r2 * np.sin(theta / 2.0), 0, r2 * np.cos(theta / 2.0)]], dtype=np.float64)
# move origin to the centre of mass
com = np.dot(masses, xyz) / np.sum(masses)
xyz = xyz - com[np.newaxis, :]
# compute inertia tensor
delta = lambda a, b: 1.0 if a==b else 0.0
imom = np.zeros((3,3), dtype=np.float64)
for i in range(3):
for j in range(3):
imom[i,j] = -np.sum(xyz[:,i] * xyz[:,j] * masses[:]) \
+ np.sum(np.sum(xyz[:,:]**2, axis=1) * masses[:]) * delta(i,j)
# diagonalize inertia tensor and compute rotational constants
diag, vec = np.linalg.eigh(imom)
A, B, C = 1.0 / diag
print("rotational constants in units 1/(au_mass * Angstrom^2):", A, B, C)
# Cartesian coordinates of atoms in the principal-axes frame
xyz = np.dot(xyz, vec)
print("Cartesian coordiates of atoms in the principal-axes frame:\n", xyz)
# convert rotational constants such that T = A*J_a^2 + B*J_b^2 + C*J_c^2
planck = 6.62606896e-27 # Plank constant in erg a second
avogno = 6.0221415e+23 # Avogadro constant
vellgt = 2.99792458e+10 # Speed of light constant in centimetres per second
convert_to_cm = planck * avogno * 1e+16 / (8.0 * np.pi * np.pi * vellgt)
print("conversion factor for rotational constants into units cm^-1:", convert_to_cm)
A, B, C = 1.0 / diag * convert_to_cm
print("rotational constants in units cm^-1:", A, B, C)
rotational constants in units 1/(au_mass * Angstrom^2): 1.5992039133513467 0.8669181785390033 0.5621696298954518
Cartesian coordiates of atoms in the principal-axes frame:
[[ 0. -0.06614561 0. ]
[ 0.75649032 0.52488941 0. ]
[-0.75649032 0.52488941 0. ]]
conversion factor for rotational constants into units cm^-1: 16.857628229754965
rotational constants in units cm^-1: 26.958785034846276 14.614184359426858 9.47684662283647
To set up the matrix representation of Hamiltonian in symmetric-top basis we need to compute the matrix elements of \(J_+^2\) and \(J_-^2\) ladder operators, as well as \(J_z^2\) and \(J^2\). However, instead of deriving and programming the explicit expressions for the matrix elements, it is sufficient to program only functions that return the result of an operator’s onto basis functions, i.e., psi2 = Jplus(psi1). The product of operators, e.g., \(J_+\cdot J_+\cdot J_+ = J_+^3\), can then be expressed as consecutive calls of operator functions, e.g., psi2 = Jplus(Jplus(Jplus(psi1))). The matrix element is then simply an overlap integral.
[14]:
# Set up functions for operators J+|psi>, J-|psi>, Jz|psi>, J^2|psi> and overlap
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)
Set up the Hamiltonian using two conventions, for near-prolate (\(z=a\)) and near-oblate (\(z=c\)) top molecules, and compute eigenvalues and eigenvectors for both cases. The resulting energies for the near-prolate and near-oblate cases must be equal, however not the eigenvectors. Print assignment of each state by the \(k_a\) and \(k_c\) quantum numbers of the \(J_z\) operator.
[15]:
# Compute rotaitonal energies and wave functions and store them in hdf5 file
import h5py
file = h5py.File("water_states.h5", "w") # open file to store the data
Jmax = 10 # maximal value of J quantum number spanned by basis
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], " %12.6f"%enr_c[istate], "ka, kc = ", abs(k_list[ind_a]), abs(k_list[ind_c]))
# store energies and wave functions in file
fgroup = file.create_group("J:"+str(J))
dset = fgroup.create_dataset('k_list', data=k_list)
dset = fgroup.create_dataset('enr', data=enr_c)
dset = fgroup.create_dataset('vec', data=vec_c)
file.close()
0 0 0.000000 0.000000 ka, kc = 0 0
1 0 24.091031 24.091031 ka, kc = 0 1
1 1 36.435632 36.435632 ka, kc = 1 1
1 2 41.572969 41.572969 ka, kc = 1 0
2 0 70.974093 70.974093 ka, kc = 0 2
2 1 79.480356 79.480356 ka, kc = 1 2
2 2 94.892369 94.892369 ka, kc = 1 1
2 3 131.926171 131.926171 ka, kc = 2 1
2 4 133.225171 133.225171 ka, kc = 2 0
3 0 138.518771 138.518771 ka, kc = 0 3
3 1 143.316754 143.316754 ka, kc = 1 3
3 2 173.927434 173.927434 ka, kc = 1 2
3 3 204.199264 204.199264 ka, kc = 2 2
3 4 210.226679 210.226679 ka, kc = 2 1
3 5 279.496300 279.496300 ka, kc = 3 1
3 6 279.709647 279.709647 ka, kc = 3 0
4 0 225.033469 225.033469 ka, kc = 0 4
4 1 227.324555 227.324555 ka, kc = 1 4
4 2 277.262887 277.262887 ka, kc = 1 3
4 3 299.536941 299.536941 ka, kc = 2 3
4 4 315.384896 315.384896 ka, kc = 2 0
4 5 377.942071 377.942071 ka, kc = 3 2
4 6 379.377116 379.377116 ka, kc = 3 1
4 7 480.549069 480.549069 ka, kc = 4 1
4 8 480.577955 480.577955 ka, kc = 4 0
5 0 330.007218 330.007218 ka, kc = 0 5
5 1 330.984120 330.984120 ka, kc = 1 5
5 2 402.751903 402.751903 ka, kc = 1 4
5 3 417.121544 417.121544 ka, kc = 2 4
5 4 448.226447 448.226447 ka, kc = 2 1
5 5 501.097085 501.097085 ka, kc = 3 3
5 6 506.385839 506.385839 ka, kc = 3 0
5 7 603.874776 603.874776 ka, kc = 4 2
5 8 604.128120 604.128120 ka, kc = 4 3
5 9 735.449590 735.449590 ka, kc = 5 1
5 10 735.453120 735.453120 ka, kc = 5 0
6 0 453.559461 453.559461 ka, kc = 0 6
6 1 453.946999 453.946999 ka, kc = 1 6
6 2 548.022973 548.022973 ka, kc = 1 5
6 3 556.066692 556.066692 ka, kc = 2 5
6 4 607.209490 607.209490 ka, kc = 0 2
6 5 648.482893 648.482893 ka, kc = 3 4
6 6 662.252959 662.252959 ka, kc = 3 1
6 7 752.550885 752.550885 ka, kc = 4 3
6 8 753.759873 753.759873 ka, kc = 4 0
6 9 883.327447 883.327447 ka, kc = 5 2
6 10 883.365500 883.365500 ka, kc = 5 3
6 11 1044.260469 1044.260469 ka, kc = 6 1
6 12 1044.260873 1044.260873 ka, kc = 6 0
7 0 595.872379 595.872379 ka, kc = 0 7
7 1 596.018910 596.018910 ka, kc = 1 7
7 2 711.567400 711.567400 ka, kc = 1 6
7 3 715.544495 715.544495 ka, kc = 2 6
7 4 790.122162 790.122162 ka, kc = 0 5
7 5 819.330163 819.330163 ka, kc = 3 5
7 6 847.406090 847.406090 ka, kc = 3 0
7 7 926.601894 926.601894 ka, kc = 4 4
7 8 930.695544 930.695544 ka, kc = 4 1
7 9 1056.648852 1056.648852 ka, kc = 5 3
7 10 1056.869847 1056.869847 ka, kc = 5 4
7 11 1216.643308 1216.643308 ka, kc = 6 2
7 12 1216.648478 1216.648478 ka, kc = 6 3
7 13 1406.989458 1406.989458 ka, kc = 7 1
7 14 1406.989503 1406.989503 ka, kc = 7 0
8 0 757.052836 757.052836 ka, kc = 0 8
8 1 757.106384 757.106384 ka, kc = 1 8
8 2 893.101745 893.101745 ka, kc = 1 7
8 3 894.895711 894.895711 ka, kc = 2 7
8 4 994.290296 994.290296 ka, kc = 0 6
8 5 1012.686628 1012.686628 ka, kc = 3 6
8 6 1060.724333 1060.724333 ka, kc = 1 1
8 7 1125.750550 1125.750550 ka, kc = 4 5
8 8 1136.544967 1136.544967 ka, kc = 4 0
8 9 1255.677608 1255.677608 ka, kc = 5 4
8 10 1256.588045 1256.588045 ka, kc = 5 5
8 11 1414.417597 1414.417597 ka, kc = 6 3
8 12 1414.452870 1414.452870 ka, kc = 6 4
8 13 1603.880403 1603.880403 ka, kc = 7 2
8 14 1603.881059 1603.881059 ka, kc = 7 3
8 15 1823.636949 1823.636949 ka, kc = 8 1
8 16 1823.636953 1823.636953 ka, kc = 8 0
9 0 937.150300 937.150300 ka, kc = 0 9
9 1 937.169376 937.169376 ka, kc = 1 9
9 2 1092.914814 1092.914814 ka, kc = 3 8
9 3 1093.673647 1093.673647 ka, kc = 2 8
9 4 1217.283817 1217.283817 ka, kc = 0 7
9 5 1227.556787 1227.556787 ka, kc = 3 7
9 6 1300.035131 1300.035131 ka, kc = 1 4
9 7 1349.384740 1349.384740 ka, kc = 4 6
9 8 1372.550710 1372.550710 ka, kc = 4 1
9 9 1480.552291 1480.552291 ka, kc = 5 3
9 10 1483.503573 1483.503573 ka, kc = 5 6
9 11 1637.876415 1637.876415 ka, kc = 6 6
9 12 1638.046289 1638.046289 ka, kc = 6 5
9 13 1826.100223 1826.100223 ka, kc = 7 5
9 14 1826.105356 1826.105356 ka, kc = 7 4
9 15 2045.043121 2045.043121 ka, kc = 8 2
9 16 2045.043200 2045.043200 ka, kc = 8 3
9 17 2294.202670 2294.202670 ka, kc = 9 1
9 18 2294.202670 2294.202670 ka, kc = 9 0
10 0 1136.185656 1136.185656 ka, kc = 0 10
10 1 1136.192318 1136.192318 ka, kc = 1 10
10 2 1311.311934 1311.311934 ka, kc = 3 9
10 3 1311.618393 1311.618393 ka, kc = 2 9
10 4 1457.869653 1457.869653 ka, kc = 0 8
10 5 1463.041567 1463.041567 ka, kc = 3 8
10 6 1562.556199 1562.556199 ka, kc = 1 5
10 7 1596.614465 1596.614465 ka, kc = 4 7
10 8 1638.539219 1638.539219 ka, kc = 2 0
10 9 1731.174563 1731.174563 ka, kc = 5 4
10 10 1739.065661 1739.065661 ka, kc = 5 7
10 11 1887.305138 1887.305138 ka, kc = 6 7
10 12 1887.949474 1887.949474 ka, kc = 6 6
10 13 2073.908402 2073.908402 ka, kc = 7 6
10 14 2073.936622 2073.936622 ka, kc = 7 5
10 15 2291.728413 2291.728413 ka, kc = 8 5
10 16 2291.729111 2291.729111 ka, kc = 8 4
10 17 2540.129357 2540.129357 ka, kc = 9 2
10 18 2540.129366 2540.129366 ka, kc = 9 3
10 19 2818.686411 2818.686411 ka, kc = 10 1
10 20 2818.686411 2818.686411 ka, kc = 10 0
We can now run same calculations using watie module
[17]:
from richmol.watie import RigidMolecule, SymtopBasis, Jxx, Jyy, Jzz
water = RigidMolecule()
water.XYZ = ("angstrom", \
"O", 0, 0, 0, \
"H", r1*np.sin(theta/2.0), 0, r1*np.cos(theta/2.0), \
"H", -r2*np.sin(theta/2.0), 0, r2*np.cos(theta/2.0) )
water.frame = "pas"
Bx, By, Bz = water.B
print("Bx, By, Bz rotational constants in cm^-1:", Bx, By, Bz)
Jmax = 10
for J in range(Jmax+1):
bas = SymtopBasis(J)
# z=c, near-oblate
H = Bx * Jxx(bas) + By * Jyy(bas) + Bz * Jzz(bas)
hmat = bas.overlap(H)
enr_c, vec = np.linalg.eigh(hmat.real)
bas_c = bas.rotate(krot=(vec.T, enr_c)) # rotate basis to the eigenvector representation
# z=a, near-prolate
H = Bz * Jxx(bas) + By * Jyy(bas) + Bx * Jzz(bas)
hmat = bas.overlap(H)
enr_a, vec = np.linalg.eigh(hmat.real)
bas_a = bas.rotate(krot=(vec.T, enr_a)) # rotate basis to the eigenvector representation
# print states energies and assignments
for istate in range(bas_c.nstates):
print( J, istate, " %12.6f"%bas_a.enr[istate], " %12.6f"%bas_c.enr[istate], \
"ka,kc =", bas_a.assign[istate][1], bas_c.assign[istate][1] )
Bx, By, Bz rotational constants in cm^-1: 26.958784982583502 14.614184327525283 9.476846602963152
0 0 0.000000 0.000000 ka,kc = 0 0
1 0 24.091031 24.091031 ka,kc = 0 1
1 1 36.435632 36.435632 ka,kc = 1 1
1 2 41.572969 41.572969 ka,kc = 1 0
2 0 70.974093 70.974093 ka,kc = 0 2
2 1 79.480356 79.480356 ka,kc = 1 2
2 2 94.892369 94.892369 ka,kc = 1 1
2 3 131.926171 131.926171 ka,kc = 2 1
2 4 133.225170 133.225170 ka,kc = 2 0
3 0 138.518771 138.518771 ka,kc = 0 3
3 1 143.316754 143.316754 ka,kc = 1 3
3 2 173.927433 173.927433 ka,kc = 1 2
3 3 204.199264 204.199264 ka,kc = 2 2
3 4 210.226679 210.226679 ka,kc = 2 1
3 5 279.496300 279.496300 ka,kc = 3 1
3 6 279.709646 279.709646 ka,kc = 3 0
4 0 225.033469 225.033469 ka,kc = 0 4
4 1 227.324555 227.324555 ka,kc = 1 4
4 2 277.262887 277.262887 ka,kc = 1 3
4 3 299.536941 299.536941 ka,kc = 2 3
4 4 315.384895 315.384895 ka,kc = 2 2
4 5 377.942071 377.942071 ka,kc = 3 2
4 6 379.377116 379.377116 ka,kc = 3 1
4 7 480.549068 480.549068 ka,kc = 4 1
4 8 480.577954 480.577954 ka,kc = 4 0
5 0 330.007218 330.007218 ka,kc = 0 5
5 1 330.984120 330.984120 ka,kc = 1 5
5 2 402.751902 402.751902 ka,kc = 1 4
5 3 417.121544 417.121544 ka,kc = 2 4
5 4 448.226446 448.226446 ka,kc = 2 1
5 5 501.097084 501.097084 ka,kc = 3 3
5 6 506.385837 506.385837 ka,kc = 3 0
5 7 603.874775 603.874775 ka,kc = 4 2
5 8 604.128119 604.128119 ka,kc = 4 3
5 9 735.449589 735.449589 ka,kc = 5 1
5 10 735.453118 735.453118 ka,kc = 5 0
6 0 453.559460 453.559460 ka,kc = 0 6
6 1 453.946998 453.946998 ka,kc = 1 6
6 2 548.022971 548.022971 ka,kc = 1 5
6 3 556.066691 556.066691 ka,kc = 2 5
6 4 607.209489 607.209489 ka,kc = 2 2
6 5 648.482892 648.482892 ka,kc = 3 4
6 6 662.252958 662.252958 ka,kc = 3 1
6 7 752.550884 752.550884 ka,kc = 4 3
6 8 753.759871 753.759871 ka,kc = 4 4
6 9 883.327445 883.327445 ka,kc = 5 2
6 10 883.365499 883.365499 ka,kc = 5 3
6 11 1044.260467 1044.260467 ka,kc = 6 1
6 12 1044.260871 1044.260871 ka,kc = 6 0
7 0 595.872378 595.872378 ka,kc = 0 7
7 1 596.018909 596.018909 ka,kc = 1 7
7 2 711.567398 711.567398 ka,kc = 1 6
7 3 715.544494 715.544494 ka,kc = 2 6
7 4 790.122161 790.122161 ka,kc = 2 5
7 5 819.330161 819.330161 ka,kc = 3 5
7 6 847.406088 847.406088 ka,kc = 3 2
7 7 926.601892 926.601892 ka,kc = 4 4
7 8 930.695542 930.695542 ka,kc = 4 1
7 9 1056.648850 1056.648850 ka,kc = 5 3
7 10 1056.869845 1056.869845 ka,kc = 5 4
7 11 1216.643305 1216.643305 ka,kc = 6 2
7 12 1216.648476 1216.648476 ka,kc = 6 3
7 13 1406.989456 1406.989456 ka,kc = 7 1
7 14 1406.989500 1406.989500 ka,kc = 7 2
8 0 757.052834 757.052834 ka,kc = 0 8
8 1 757.106382 757.106382 ka,kc = 1 8
8 2 893.101743 893.101743 ka,kc = 1 7
8 3 894.895709 894.895709 ka,kc = 2 7
8 4 994.290294 994.290294 ka,kc = 0 6
8 5 1012.686626 1012.686626 ka,kc = 3 6
8 6 1060.724331 1060.724331 ka,kc = 1 1
8 7 1125.750548 1125.750548 ka,kc = 4 5
8 8 1136.544965 1136.544965 ka,kc = 4 6
8 9 1255.677606 1255.677606 ka,kc = 5 4
8 10 1256.588042 1256.588042 ka,kc = 5 5
8 11 1414.417594 1414.417594 ka,kc = 6 3
8 12 1414.452867 1414.452867 ka,kc = 6 4
8 13 1603.880400 1603.880400 ka,kc = 7 2
8 14 1603.881056 1603.881056 ka,kc = 7 3
8 15 1823.636945 1823.636945 ka,kc = 8 1
8 16 1823.636950 1823.636950 ka,kc = 8 2
9 0 937.150298 937.150298 ka,kc = 0 9
9 1 937.169374 937.169374 ka,kc = 1 9
9 2 1092.914812 1092.914812 ka,kc = 3 8
9 3 1093.673645 1093.673645 ka,kc = 2 8
9 4 1217.283814 1217.283814 ka,kc = 0 7
9 5 1227.556785 1227.556785 ka,kc = 3 7
9 6 1300.035128 1300.035128 ka,kc = 1 4
9 7 1349.384737 1349.384737 ka,kc = 4 6
9 8 1372.550707 1372.550707 ka,kc = 4 1
9 9 1480.552288 1480.552288 ka,kc = 5 3
9 10 1483.503570 1483.503570 ka,kc = 5 6
9 11 1637.876412 1637.876412 ka,kc = 6 6
9 12 1638.046286 1638.046286 ka,kc = 6 5
9 13 1826.100219 1826.100219 ka,kc = 7 5
9 14 1826.105353 1826.105353 ka,kc = 7 4
9 15 2045.043117 2045.043117 ka,kc = 8 2
9 16 2045.043196 2045.043196 ka,kc = 8 3
9 17 2294.202665 2294.202665 ka,kc = 9 1
9 18 2294.202666 2294.202666 ka,kc = 9 2
10 0 1136.185654 1136.185654 ka,kc = 0 10
10 1 1136.192315 1136.192315 ka,kc = 1 10
10 2 1311.311931 1311.311931 ka,kc = 3 9
10 3 1311.618390 1311.618390 ka,kc = 2 9
10 4 1457.869650 1457.869650 ka,kc = 4 8
10 5 1463.041564 1463.041564 ka,kc = 3 8
10 6 1562.556196 1562.556196 ka,kc = 1 5
10 7 1596.614462 1596.614462 ka,kc = 4 7
10 8 1638.539216 1638.539216 ka,kc = 2 2
10 9 1731.174560 1731.174560 ka,kc = 5 4
10 10 1739.065658 1739.065658 ka,kc = 5 7
10 11 1887.305134 1887.305134 ka,kc = 6 7
10 12 1887.949470 1887.949470 ka,kc = 6 6
10 13 2073.908398 2073.908398 ka,kc = 7 6
10 14 2073.936618 2073.936618 ka,kc = 7 5
10 15 2291.728408 2291.728408 ka,kc = 8 5
10 16 2291.729106 2291.729106 ka,kc = 8 4
10 17 2540.129352 2540.129352 ka,kc = 9 2
10 18 2540.129361 2540.129361 ka,kc = 9 3
10 19 2818.686405 2818.686405 ka,kc = 10 1
10 20 2818.686405 2818.686405 ka,kc = 10 2
Remark: for some reason, rotational constants calculated explicitly above and obtained from watie have slightly differrent values (need to find the reason), which is reflected in the assignment of rotaitonal states for high \(J\).
[ ]: