Skip to content
scipython.com

E6.21Modelling the distribution of \(\mathrm{^{13}C}\) atoms in \(\mathrm{C_{60}}\)

There are two stable isotopes of carbon, \(\mathrm{^{12}C}\) and \(\mathrm{^{13}C}\) (the radioactive \(\mathrm{^{14}C}\) nucleus is present in nature in only trace amounts of the order of parts per trillion). Taking the abundance of \(\mathrm{^{13}C}\) to be \(x=0.0107\) (i.e. about 1%), we will calculate the relative amounts of buckminster fullerene, \(\mathrm{C_{60}}\), with exactly zero, one, two, three and four \(\mathrm{^{13}C}\) atoms. (This is important in nuclear magnetic resonance studies of fullerenes, for example, because only the \(\mathrm{^{13}C}\) nucleus is magnetic and so detectable by NMR).

The number of \(\mathrm{^{13}C}\) atoms in a population of carbon atoms sampled at random from a population with natural isotopic abundance follows a binomial distribution: the probability that, out of \(n\) atoms, \(m\) will be \(\mathrm{^{13}C}\) (and therefore \(n-m\) will be \(\mathrm{^{12}C}\)) is

\[p_m(n) = \binom{n}{m} x^m (1-x)^{n-m}.\]

We can, of course, calculate \(p_m(60)\) exactly from this formula for \(0 \le m \le 4\), but we can also simulate the sampling with the np.random.binomial method:

c13inc60
import numpy as np

n, x = 60, 0.0107
mmax = 4
m = np.arange(mmax + 1)

# Estimate the abundances by random sampling from the binomial distribution
ntrials = 10000
pbin = np.empty(mmax + 1)
for r in m:
    pbin[r] = np.sum(np.random.binomial(n, x, ntrials) == r) / ntrials

# Calculate and store the binomial coefficients nCm
nCm = np.empty(mmax + 1)
nCm[0] = 1
for r in m[1:]:
    nCm[r] = nCm[r - 1] * (n - r + 1) / r
# The "exact" answer from binomial distribution
p = nCm * x**m * (1 - x) ** (n - m)

print("Abundances of C60 as (13C)[m](12C)[60-m]")
print('m   "Exact"    Estimated')
print("-" * 24)
for r in m:
    print(f"{r:1d}   {p[r]:6.4f}      {pbin[r]:6.4f}")

For each value of r in the array m, we sample a large number of times (ntrial) from the binomial distribution described by \(n=60\) and probability, \(x=0.0107\). The comparison of these sample values with a given value of r yields a boolean array which can be summed (remembering that True evaluates to 1 and False evaluates to 0); division by ntrials then gives an estimate of the probability of exactly r atoms being of type \(\mathrm{^{13}C}\) and the remainder of type \(\mathrm{^{12}C}\).

The explicit loop over m could be removed by creating an array of shape (ntrials, mmax+1) containing all the samples, and summing over the first axis of this array in the comparison with the m array:

samples = np.random.binomial(n, x, (ntrials, mmax+1))
pbin = np.sum(samples == m, axis=0) / ntrials

The abundances of \(\mathrm{^{13}C_m^{12}C_{60-m}}\) produced by our program are given as the following output.

Abundances of C60 as (13C)[m](12C)[60-m]
m   "Exact"    Estimated
------------------------
0   0.5244      0.5199
1   0.3403      0.3348
2   0.1086      0.1093
3   0.0227      0.0231
4   0.0035      0.0031

That is, almost 48% of \(\mathrm{C_{60}}\) molecules contain at least one magnetic nucleus.