Modelling 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:

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('{:1d}   {:6.4f}      {:6.4f}'.format(r, p[r], pbin[r]))

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.