Isotopologues

In the previous section we successfully assigned and modeled the rotational spectrum of the main isotopologue of OCS. Now we want to see which other isotopologues we can find in the spectrum. Isotopologues are very interesting as they are excellent tracers of isotope ratios (e.g., for astronomers) and can be used to determine the structure of the molecule. We will start with a small little exercise. Research the natural abundancees and masses of the different oxygen, carbon, and sulfur isotopes. Then calculate the natural abundancees of the resulting isotopologues and estimate their rotational constants $B$ in the same way we did for the main isotopologue previously.

Solution

We copy the values for the $r_{CO}$ and $r_{CS}$ distances from the previous exercise and research the masses and abundancees of the most prevalent isotopes. The abundancees of the isotopologues are then calculated as the product of the abundancees of their isotopes. For the $B$ rotational constant, we use the same equation as in the previous exercise but with the masses for the isotopologues. See the following code for a possible Python implementation.

			
import numpy as np 

r_CO = 1.1308051672090703e-10
r_CS = 1.5375286554975936e-10
m0 = 1.66054e-27
h = 6.626e-34

O_isotopes = [[16, 0.998,   15.994914619257],
				[17, 0.00038, 16.999131755953],
				[18, 0.00205, 17.999159612136]]

C_isotopes = [[12, 0.989, 12],
				[13, 0.0106, 13.003354835336]]

S_isotopes = [[32, 0.948, 31.9720711744],
				[33, 0.0076, 32.9714589099],
				[34, 0.0437, 33.96786701],
				[36, 0.0002, 35.96708070]]


Bs = []
for nO, pO, mO in O_isotopes:
	for nC, pC, mC in C_isotopes:
			for nS, pS, mS in S_isotopes:
				masses = np.array((mO, mC, mS))
				positions = np.array((0, r_CO, r_CO + r_CS))                
				c = np.sum(masses * positions) / np.sum(masses)
				positions -= c
				
				I = np.sum(masses * m0 * positions**2)
				B = h / (8 * np.pi**2 * I)
				abundance = pO * pC * pS * 100
				Bs.append(((nO, nC, nS), B/1E6, abundance))


Bs = sorted(Bs, key=lambda x: x[-1], reverse=True)
output = ['| mO | mC | mS | B [MHz] | Abundance|', '|----|----|----|---------|----------|']
for (nO, nC, nS), B, percentage in Bs:
	output.append(f'| {nO} | {nC} | {nS} | {B:.2f} | {percentage:7.4f}% |')

print("\n".join(output))

The results are summarized in the following table:

mO	mC	mS	B [MHz]	Abundance
16	12	32	6322.22	93.5697%
16	12	34	6167.36	4.3133%
16	13	32	6301.34	1.0029%
16	12	33	6242.45	0.7501%
18	12	32	5930.97	0.1922%
16	13	34	6144.89	0.0462%
17	12	32	6116.73	0.0356%
16	12	36	6028.68	0.0197%
18	12	34	5780.04	0.0089%
16	13	33	6220.75	0.0080%
18	13	32	5916.16	0.0021%
17	12	34	5963.94	0.0016%
18	12	33	5853.22	0.0015%
17	13	32	6099.16	0.0004%
17	12	33	6038.03	0.0003%
16	13	36	6004.73	0.0002%
18	13	34	5763.87	0.0001%
18	12	36	5644.86	0.0000%
17	13	34	5944.91	0.0000%
18	13	33	5837.72	0.0000%
17	12	36	5827.10	0.0000%
17	13	33	6019.71	0.0000%
18	13	36	5627.42	0.0000%
17	13	36	5806.70	0.0000%

Obviously we start with the most abundant isotopologue (OC$^{34}$S at 4%). Its estimated rotational constant is $B=6167.36 \text{ MHz}$ but we have some additional information. From our analysis of the main isotopologue, we know that the estimated rotational constant was a little bit too large. As the rotational constant depends only on the geometry and the masses, probably our geometry was a little bit off (which is not very surprising given the estimations we made). We can calculate a correction factor from the experimental value and the initial estimated value for $B$. This should yield $B_\text{exp} / B_\text{est} \approx 0.962$. Therefore we should multiply all the estimated $B$ values with this correction factor.

Solution

mO	mC	mS	B [MHz]	Abundance
16	12	32	6081.98	93.5697%
16	12	34	5933.00	4.3133%
16	13	32	6061.89	1.0029%
16	12	33	6005.24	0.7501%
18	12	32	5705.59	0.1922%
16	13	34	5911.38	0.0462%
17	12	32	5884.29	0.0356%
16	12	36	5799.59	0.0197%
18	12	34	5560.40	0.0089%
16	13	33	5984.37	0.0080%
18	13	32	5691.34	0.0021%
17	12	34	5737.31	0.0016%
18	12	33	5630.80	0.0015%
17	13	32	5867.40	0.0004%
17	12	33	5808.58	0.0003%
16	13	36	5776.55	0.0002%
18	13	34	5544.84	0.0001%
18	12	36	5430.35	0.0000%
17	13	34	5719.00	0.0000%
18	13	33	5615.89	0.0000%
17	12	36	5605.67	0.0000%
17	13	33	5790.96	0.0000%
18	13	36	5413.58	0.0000%
17	13	36	5586.05	0.0000%

Next, we create a model for the OC$^{34}$S isotopologue. For this we can copy the main isotopologue *.par file and only change the $B$ value as it is often better to fix the higher-order constants to the main isotopologue value than setting them to zero. You can find the respective *.var file here (in case you want to compare your version). We will use the same *.int file for now and neglect any small changes in the dipole moment. Create predictions for the transitions by running SPCAT and load them together with the spectrum into LLWP. You should immediately see a pattern of experimental lines close to the predicted transitions.

The model for the isotopologue scaled by the main isotopologue ratio of the experimental and estimated $B$ value reproduces the spectrum of OC$^{34}$S already surprisingly well. The deviations are easy to follow in the Loomis-Wood plot which makes the assignments efficient and reliable.

As an exercise, create a second set of predictions without the higher-order parameters and compare the two Loomis-Wood plots with each other. You should see that fixing the constants to the main isotopologue value results in more accurate predictions.

Same Loomis-Wood plot as for the previous figure but without the higher-order constants (fixed to the main isotopologue values). The deviations are much more pronounced here (note the different widths of the two plots: 20 MHz for the previous one and 80 MHz for this one).

Assign-Fit-Predict

Now repeat the same steps as for the main isotopologue:

Assign the correct quantum numbers to the experimental lines in LLWP
Save the assignments to a *.lin file
Create a *.par file from the *.var file (update the uncertainties, number of lines)
Run SPFIT
Check the relative uncertainties of the parameters (are all parameters defined, if not, remove them, repeat)
Run SPCAT
Check the residuals

Similarly to the main isotopologue you can determine sensible uncertainties for your new assignments and change the uncertainties in the *.lin file accordingly. Then you can search for OC$^{34}$S in the CDMS and add literature data to your model. Check if you can/should add any higher-order parameters now. Finally, compare your model to the model found in the CDMS.

No Model - No Problem

After having successfully modeled OCS and OC$^{34}$S, we continue our search of more isotopologues in the spectrum. However, this time we will search for patterns in the Loomis-Wood plots and forego initial predictions.

Load the OCS spectrum and the *.cat and *.int files of the two isotopologues that we have modelled already. Loading the *.lin files marks all already assigned lines with a blue star. This makes it easy to exclude any lines we have already assigned. Then we can use the predictions for the main isotopologue and search for patterns around it. However, there is one problem. We have loaded two *.cat and *.int files with the same quantum numbers. LLWP will use the first matching prediction for the reference series (which are the OC$^{34}$S in this case). To limit the selection to a single file, go to the Reference Series window and click the All button behind the File: label. In the dialog that opens up, choose the file with the main isotopologue data to use as our reference file. This is very useful when working with many single state fits (as we are doing here).

Loomis-Wood plot after loading the *.cat and *.lin files of the two already analyzed isotopologues. The red sticks indicate our reference series (here the main isotopologue data) and purple indicates any other predicted transitions. The blue stars indicate an assignment.

The resulting Loomis-Wood plot makes it very tough to see any additional isotopologues as each row is scaled to its maximum intensity. Thus any less intense lines are hard to spot. We can change this by opening the Scaling Window via the menu (Plot > Scaling Window). Then choose in there the Custom option for the y-axis. New input fields will appear which specify the min and max values of the y-axis (exclusive of a user-defined margin) and the scaling factor between the predicted and experimental intensities. Choose 6.5e2 for the scaling and 0 and 0.001 for the min and max values, respectively. Now two additional series should be easy to see in the Loomis-Wood plot.

Two so-far unassigned series are visible in the Loomis-Wood plot when adjusting the scale.

Now assign one of the two series, save the assignments to a new *.lin file and fit a model to it. Once your models for the two series are finalized, compare the resulting $B$ values with the estimated $B$ values. Which two isotopologues did you assign?

Solution

You have assigned OC$^{33}$S and O$^{13}$CS. OC$^{33}$S is already evident from the spectrum as its series is right between OCS and OC$^{34}$S. The transitions of O$^{13}$CS (and hence also its rotational constant $B$) are very close to those of the main isotopologue, as the carbon atom is very close to the center of mass.

This little exercise is meant to show that we can find similar fingerprints (of isotopologues or vibrationally excited states) by searching for similar patterns close-by in the Loomis-Wood plot. Conversely, if you spot any of these patterns in Loomis-Wood plots, you know now what most likely their origin is.

Linear Molecules in SPFIT/SPCAT

Vibrational States