Molecule input
One can define molecule using one of the following ways
Provide Cartesian coordinates of atoms, dipole moment, polarizability, etc., referring to the same (user-defined) coordinate frame. These can be almost routinely obtained from quantum chemical calculations.
Provide spectroscopic constants (\(A\), \(B\), \(C\), \(\Delta_{J}\), \(\Delta_{J,K}\) …) and dipole moment, polarizability, etc., in the coordinate frame of principal axes of inertia (PAI).
Provide spectroscopic constants (\(A\), \(B\), \(C\), \(\Delta_{J}\), \(\Delta_{J,K}\) …) and dipole moment, polarizability, etc., which are not in the PAI frame. In this case, Cartesian coordinates of atoms defined with respect to the same frame as dipole moment, polarizability, etc., must be provided, although they will be used only for ‘rotating’ the coordinate frame to the PAI system.
[37]:
from richmol.rot import Molecule, mol_tensor
Cartesian data
Here is an example for water molecule using data obtained from a quantum chemical calculation
[38]:
water = Molecule()
# Cartesian coordinates of atoms
water.XYZ = ("bohr",
"O", 0.00000000, 0.00000000, 0.12395915,
"H", 0.00000000, -1.43102686, -0.98366080,
"H", 0.00000000, 1.43102686, -0.98366080)
# dipole moment (au)
water.dip = [0, 0, -0.7288]
# polarizability tensor (au)
water.pol = [[9.1369, 0, 0], [0, 9.8701, 0], [0, 0, 9.4486]]
print("Masses and Cartesian coordinates of atoms")
for atom in water.XYZ:
print(atom['label'], atom['mass'], atom['xyz'])
Masses and Cartesian coordinates of atoms
O 15.99491462 [0. 0. 0.06559636]
H 1.0078250322 [ 0. -0.7572668 -0.52053088]
H 1.0078250322 [ 0. 0.7572668 -0.52053088]
By default, calculations are carried out for the main isotopologue. To specify a non-standard isotope, put the corresponding isotope number next to the atom label, for example, “O18” for oxygen-18 or “H2” for deuterium
[39]:
# example of D2^{18}O
D2_18O = Molecule()
D2_18O.XYZ = ("bohr",
"O18", 0.00000000, 0.00000000, 0.12395915,
"H2", 0.00000000, -1.43102686, -0.98366080,
"H2", 0.00000000, 1.43102686, -0.98366080)
print("Masses and Cartesian coordinates of atoms")
for atom in D2_18O.XYZ:
print(atom['label'], atom['mass'], atom['xyz'])
Masses and Cartesian coordinates of atoms
O18 17.999159613 [0. 0. 0.06559636]
H2 2.0141017781 [ 0. -0.7572668 -0.52053088]
H2 2.0141017781 [ 0. 0.7572668 -0.52053088]
One can also read/store the geometry from/to XYZ file
[40]:
# store Cartesian coordinates of atoms into XYZ file
water.store_xyz("water.xyz", "some comment line")
# read Cartesian coordinates of atoms from XYZ file
water2 = Molecule()
water2.XYZ = "water.xyz"
print("Masses and Cartesian coordinates of atoms")
for atom in water2.XYZ:
print(atom['label'], atom['mass'], atom['xyz'])
Masses and Cartesian coordinates of atoms
O 15.99491462 [0. 0. 0.06559636]
H 1.0078250322 [ 0. -0.7572668 -0.52053088]
H 1.0078250322 [ 0. 0.7572668 -0.52053088]
Coordinate frame
The molecular frame embedding, i.e., orientation of the \(x\), \(y\), \(z\) axes with respect to molecule, can be altered using frame. The tensor properties dip for dipole moment and pol for polarizability tensor (and few others) will be automatically rotated to a new coordinate frame whenever the latter is changed.
The frame can be defined by its name, e.g., frame="ipas". One can use the name of a Molecule attribute containing a (3x3) arbitrary rotation matrix, e.g., water.rotmat=[[0,0,1],[1,0,0],[0,1,0]]; frame="rotmat". A third option is to combine a linear operation with the name of an attribute, e.g., water.pol = [[9.1369,0.3,-0.01],[0.3, 9.8701,0],[-0.01,0,9.4486]]; frame="diag(pol)" will rotate frame to the principal axes of polarizability (i.e., frame where water.pol
tensor becomes diagonal). Axes permutation is also a frame rotation, e.g., frame="zxy" performs (123) permutation of the axes. Several examples are shown below.
[41]:
# change frame to inertial principal axes system (ipas)
water.frame = "ipas"
# or equivalently
water.frame = "diag(inertia)" # reads as: frame where the inertia tensor is diagonal
print("coordinates of atoms\n", water.XYZ)
print("dipole moment\n", water.dip)
print("polarizability\n", water.pol)
print("inertia tensor\n", water.inertia)
coordinates of atoms
[('O', 15.99491462, [ 0. , 0.06559636, 0. ])
('H', 1.00782503, [-0.7572668 , -0.52053088, 0. ])
('H', 1.00782503, [ 0.7572668 , -0.52053088, 0. ])]
dipole moment
[ 0. -0.7288 0. ]
polarizability
[[9.8701 0. 0. ]
[0. 9.4486 0. ]
[0. 0. 9.1369]]
inertia tensor
[[ 0.61496945 -0. -0. ]
[-0. 1.1558806 -0. ]
[-0. -0. 1.77085004]]
Multiple frame operations can be combined together, for example, we can rotate to the inertial principal axes system and then permute \(x\) and \(z\) axes
[42]:
water.frame = "ipas"
water.frame = "zyx"
# or equivalently in one line (with operations performed form right to left)
water.frame = "zyx,ipas"
print("coordinates of atoms\n", water.XYZ)
print("dipole moment\n", water.dip)
print("polarizability\n", water.pol)
print("inertia tensor\n", water.inertia)
coordinates of atoms
[('O', 15.99491462, [ 0. , 0.06559636, 0. ])
('H', 1.00782503, [ 0. , -0.52053088, -0.7572668 ])
('H', 1.00782503, [ 0. , -0.52053088, 0.7572668 ])]
dipole moment
[ 0. -0.7288 0. ]
polarizability
[[9.1369 0. 0. ]
[0. 9.4486 0. ]
[0. 0. 9.8701]]
inertia tensor
[[ 1.77085004 -0. -0. ]
[-0. 1.1558806 -0. ]
[-0. -0. 0.61496945]]
The principal axes system can be defined with respect to any rank-2 symmetric matrix. For example, in many cases it is convenient to choose molecular frame such that polarizability tensor becomes diagonal
[43]:
water.frame = "diag(pol)" # reads as: frame where tensor 'pol' is diagonal
print("coordinates of atoms\n", water.XYZ)
print("dipole moment\n", water.dip)
print("polarizability\n", water.pol)
coordinates of atoms
[('O', 15.99491462, [ 0. , 0.06559636, 0. ])
('H', 1.00782503, [ 0. , -0.52053088, -0.7572668 ])
('H', 1.00782503, [ 0. , -0.52053088, 0.7572668 ])]
dipole moment
[ 0. -0.7288 0. ]
polarizability
[[9.1369 0. 0. ]
[0. 9.4486 0. ]
[0. 0. 9.8701]]
We can also define a custom rotation matrix or principal axis matrix, as demonstrated below
[44]:
import numpy as np
import scipy.linalg as la
# generate random matrix
mat = np.random.rand(3,3)
# random rotation matrix
water.custom_rot = la.expm((mat - mat.T)/2)
# random symmetric matrix
water.custom_pam = (mat + mat.T)/2
# use 'custom_rot' (orthogonal) matrix to rotate the frame
water.frame = "custom_rot"
# alternatively use 'custom_pam' (symmetric) matrix as principal axes matrix
water.frame = "diag(custom_pam)"
Multiple occurences of frame assignments lead to the accumulation of corresponding frame transformations. Use frame=None or frame="None" to reset the frame to the one defined by the input Cartesian coordinates of atoms
[45]:
water.frame = "diag(inertia)" # rotate to ipas
water.frame = "zxy" # then permute axes (123)
# now we want to permute axes (23) in the original frame, reset frame rotation
water.frame = None
water.frame = "xzy"
# or equivalently in one line
water.frame = "xzy,None"
print(water.XYZ)
[('O', 15.99491462, [ 0. , 0.06559636, 0. ])
('H', 1.00782503, [ 0. , -0.52053088, -0.7572668 ])
('H', 1.00782503, [ 0. , -0.52053088, 0.7572668 ])]
Only certain molecule properties (such as dip, pol, inertia, XYZ, and few others) are automatically rotated when the frame is altered. This is demonstrated in the following example
[46]:
water.vec = [1,2,3]
water.mat = [[1,2,3],[4,5,6],[7,8,9]]
water.frame = "zyx,null"
# axes permutation (13) does not affect vec and mat, but dip
print("vec\n", water.vec)
print("mat\n", water.mat)
print("dip\n", water.dip)
vec
[3. 2. 1.]
mat
[[9. 8. 7.]
[6. 5. 4.]
[3. 2. 1.]]
dip
[-0.7288 0. 0. ]
To declare a custom Cartesian tensor which is dynamically rotated to a new frame, use mol_tensor function
[47]:
water.vec = mol_tensor([1,2,3])
water.mat = mol_tensor([[1,2,3],[4,5,6],[7,8,9]])
water.frame = "zyx,null"
# axes permutation (13) now affects also vec and mat
print("vec\n", water.vec)
print("mat\n", water.mat)
print("dip\n", water.dip)
# we can then change vec and mat values, they will still be dynamically rotated
water.vec = [7,8,9]
water.mat = [[11,12,13],[14,15,16],[17,18,19]]
print("\nnew vec\n", water.vec)
print("new mat\n", water.mat)
vec
[3. 2. 1.]
mat
[[9. 8. 7.]
[6. 5. 4.]
[3. 2. 1.]]
dip
[-0.7288 0. 0. ]
new vec
[9. 8. 7.]
new mat
[[19. 18. 17.]
[16. 15. 14.]
[13. 12. 11.]]
Rotational constants
The rotational constants \(B_x\), \(B_y\), \(B_z\) (in units of cm\(^{-1}\)) can be calculated from the input molecular geometry using B. This will work only when the molecular frame is set to the inertial principal axes system, i.e., when the inertia tensor is diagonal
[48]:
water.frame = "diag(pol)" # frame where polarizability tensor is diagonal
print(water.B) # works, because inertia tensor happens to be diagonal in the selected frame
water.frame = "diag(inertia)"
print(water.B) # the rotational constants in the inertia frame will be sorted in ascending order
[9.519512542010801, 14.584230612810305, 27.412141069544163]
[27.412141069544163, 14.584230612810305, 9.519512542010801]
B can be also used to assign user-defined values to the rotational constants, for example, experimental values. Once assigned, B will always return the user-defined values. To access the values of rotational constants calculated from molecular geometry, use B_geom
[49]:
print("water.B ", water.B)
water.B = (27.877, 14.512, 9.285) # user-defined rotational constants in units of cm^-1
print("water.B ", water.B)
print("water.B_geom ", water.B_geom)
water.B [27.412141069544163, 14.584230612810305, 9.519512542010801]
water.B [27.877, 14.512, 9.285]
water.B_geom [27.412141069544163, 14.584230612810305, 9.519512542010801]
Important: If rotational constants were assigned to user-defined values, the rotational Hamiltonian is built using these values and not the molecular geometry input
For linear and spherical-top molecules, use B to print and assign the rotational constant
[50]:
ocs = Molecule()
ocs.B = 0.20286 # experimental value
print(ocs.B)
ocs.XYZ = ("angstrom", "C", 0, 0, 0, "S", 0, 0, -1.56, "O", 0, 0, 1.16)
print(ocs.B_geom) # calculated value
print("molecule is linear ?:", ocs.linear()) # True if molecule is linear
[0.20286, 0.20286, 0]
[0.20317859211971495, 0.20317859211971495, 0]
molecule is linear ?: True
If molecular geometry is specified, the user-assigned rotational constants are always verified to agree within the 5% difference with the values calculated from the input molecular geometry. Here is an example of an error when the assigned and calculated constants do not match
[51]:
water2 = Molecule()
water2.B = (27.877, 9.285, 14.512) # user-defined rotational constants in units of cm^-1
water2.XYZ = ("bohr",
"O", 0.00000000, 0.00000000, 0.12395915,
"H", 0.00000000, -1.43102686, -0.98366080,
"H", 0.00000000, 1.43102686, -0.98366080)
water2.frame = 'ipas' # inertia principal axes system
print("inertia frame")
print("calculated Bx, By, Bz", water2.B_geom)
print("input Bx, By, Bz ", water2.B)
inertia frame
calculated Bx, By, Bz [27.412141069544163, 14.584230612810305, 9.519512542010801]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-51-a0f2da6fa694> in <module>
11 print("inertia frame")
12 print("calculated Bx, By, Bz", water2.B_geom)
---> 13 print("input Bx, By, Bz ", water2.B)
~/.local/lib/python3.8/site-packages/richmol-0.1a1-py3.8-linux-x86_64.egg/richmol/rot/molecule.py in B(self)
316 """
317 try:
--> 318 self.check_B()
319 return self.B_exp
320 except AttributeError:
~/.local/lib/python3.8/site-packages/richmol-0.1a1-py3.8-linux-x86_64.egg/richmol/rot/molecule.py in check_B(self)
346 err = "\n".join( abc + " %12.6f"%e + " %12.6f"%c + " %12.6f"%(e-c) \
347 for abc,e,c in zip(("Bx","By","Bz"), B_exp, B_geom) )
--> 348 raise ValueError(f"input rotational constants disagree much with geometry\n" + \
349 f" inp calc inp-calc\n" + err) from None
350 return
ValueError: input rotational constants disagree much with geometry
inp calc inp-calc
Bx 27.877000 27.412141 0.464859
By 9.285000 14.584231 -5.299231
Bz 14.512000 9.519513 4.992487
This can be fixed, for example, if we swap \(y\) and \(z\) axes
[52]:
water2.frame = 'xzy,ipas'
print("inertia frame")
print("calculated Bx, By, Bz", water2.B_geom)
print("input Bx, By, Bz ", water2.B)
inertia frame
calculated Bx, By, Bz [27.412141069544163, 9.519512542010801, 14.584230612810305]
input Bx, By, Bz [27.877, 9.285, 14.512]
Centrifugal distortion constants
The following Watson-type asymmetric top Hamiltonians in the \(A\) or \(S\) standard reduced form are implemented (they can also be easily modified and extended to incorporate higher-order terms)
A-form: \(H_A = H_\text{rigrot} - \Delta_{J} J^{4} - \Delta_{JK} J^{2} J_{z}^{2} - \Delta_{K} J_{z}^{4} - \frac{1}{2} [ \delta_{J} J^{2} + \delta_{K} J_{z}^{2}, J_{+}^{2} + J_{-}^{2} ]_{+} + H_{J} J^{6} + H_{JK} J^{4} J_{z}^{2} + H_{KJ} J^{2} J_{z}^{4} + H_{K} J_{z}^{6} + \frac{1}{2} [ \phi_{J} J^{4} + \phi_{JK} J^{2} J_{z}^{2} + \phi_{K} J_{z}^{4}, J_{+}^{2} + J_{-}^{2} ]_{+}\)
S-form: \(H_S = H_\text{rigrot} - \Delta_{J} J^{4} - \Delta_{JK} J^{2} J_{z}^{2} - \Delta_{K} J_{z}^{4} + d_{1} J^{2} (J_{+}^{2} + J_{-}^{2}) + d_{2} (J_{+}^{4} + J_{-}^{4}) + H_{J} J^{6} + H_{JK} J^{4} J_{z}^{2} + H_{KJ} J^{2} J_{z}^{4} + H_{K} J_{z}^{6} + h_{1} J^{4} (J_{+}^{2} + J_{-}^{2}) + h_{2} J^{2} (J_{+}^{4} + J_{-}^{4}) + h_{3} (J_{+}^{6} + J_{-}^{6})\)
The type of Watson reduction can be specified using watson. The centrifugal constants are defined as molecule attributes DeltaJ, DeltaJK, DeltaK, deltaJ, deltaK, HJ, HJK, HKJ, HK, phiJ, phiJK, phiK, d1, d2, h1, h2, h3. Here is an example of setting Watson A-type Hamiltonian for water using values from NIST
[36]:
from richmol.rot.molecule import Molecule
from richmol.rot.solution import solve
from richmol.convert_units import MHz_to_invcm
import sys
# NIST values of rotational constants for H2O (in MHz)
water = Molecule()
water.B = MHz_to_invcm(835840.288, 435351.717, 278138.700)
water.watson = 'watson_a' # or 'watson_s' for S-type
water.DeltaJ = 37.59422 * MHz_to_invcm()
water.DeltaJK = -172.9128 * MHz_to_invcm()
water.DeltaK = 973.29052 * MHz_to_invcm()
water.deltaJ = 15.210402 * MHz_to_invcm()
water.deltaK = 41.0502 * MHz_to_invcm()
water.HJ = 1.56556e-2 * MHz_to_invcm()
water.HJK = -4.2081e-2 * MHz_to_invcm()
water.HKJ = -5.09508e-1 * MHz_to_invcm()
water.HK = 3.733028 * MHz_to_invcm()
water.phiJ = 7.79579e-3 * MHz_to_invcm()
water.phiJK = -2.5165e-2 * MHz_to_invcm()
water.phiK = 1.0971 * MHz_to_invcm()
[ ]: