Generate raw combinations of elements or species

This notebook aims to demonstrate the sheer scale of the combinatorial explosion. Based on some assumtions about how many elements (103) or species - elements in oxidation states - (403) there are, it calculates how many binary, ternary, quaternary... combinations there can be, within different integers stoichiometry limits.

this notebook contains ZERO chemistry.

The numbers that are produced match roughly to the raw combinations (before chemical filters) numbers obtained in Table 1 of the publication Computational Screeningf of All Stoichiometric Inorganic Materials30155-3) (D. W. Davies et. al, Chem, 2016). Small deviations from the numbers are due to changes in what we class as an "accessible" oxidation state for some elements.

N.B. It can take a minute or two to calculate quaternaries on a desktop/laptop computer.


In [1]:
# imports 
import itertools
from fractions import gcd

In [2]:
# Functions
def sum_iter(iterator):
    return sum(1 for _ in iterator)


def _gcd_recursive(*args):
    """
    Get the greatest common denominator among any number of ints
    """
    if len(args) == 2:
        return gcd(*args)
    else:
        return gcd(args[0], _gcd_recursive(*args[1:]))


def _is_irreducible(stoichs):
    for i in range(1, len(stoichs)):
        if gcd(stoichs[i-1], stoichs[i]) == 1:
            return True
    else:
        return False

# MAIN
# >>>     Set to 103 for combinations of elements            <<<
# >>> or 403 for combinations of element in oxidation states <<<
search_space = 403


def main():

    compound_names = {2:'binary', 3: 'ternary',
                      4: 'quaternary', 5:'quinternary'}

    print ("In a search space of {0} elements:".format(search_space))
    print ("")

    def count_combinations(n):
        c = sum_iter(itertools.combinations(range(search_space), n))
        print ("Number of unique combinations of {0} elements: {1}".format(n, c))
        return c

    combinations = {n: count_combinations(n) for n in range(2,5)}


    print("")

    # Count all the combinations of inequivalent ratios
    # e.g. for ternaries up to x=2: (1,1,1), (1,1,2), (1,2,1),
    #                               (1,2,2), (2,1,1), (2,1,2), (2,1,2)
    def ineq_ratios_with_coeff(n_species, max_coeff):
        print ("Number of inequivalent {0} ratios with max coefficient {1}: ".format(
                   compound_names[n_species], max_coeff)),
        n_ratios = sum_iter(filter(
            _is_irreducible,   # Filter operation cuts equivalent
                               # stoichs e.g. accept (1,1), reject (2,2)
            itertools.product(*( # All stoich combinations, e.g.
                                  # Product of (1,2), (1,2) is
                                  # (1,1), (2,1), (2,1), (2,2)
                                  #
                n_species * (tuple(range(1,max_coeff+1)),)
                                    # Multiply nested tuples to form
                                    # initial set of possibilities
                                    # (1,2,..), (1,2,..),...
                                ))))

        print ("{0}".format(n_ratios))
        return n_ratios

    ratios = {n: {x: ineq_ratios_with_coeff(n, x) for x in range(2,9,2)} for n in range(2,5)}
    print ("")


    def unique_compounds(n_species, max_coeff):
        print ("Unique {0} compounds (no stability constraints) with max coefficient {1}: ".format(
        compound_names[n_species], max_coeff)),
        n_compounds = combinations[n_species] * ratios[n_species][max_coeff]
        print ("{0:5.3e}".format(n_compounds))
        return n_compounds

    compounds = {n: {x: unique_compounds(n, x) for x in range(4,9,2)} for n in range(2,5)}

if __name__ == '__main__':
    main()


In a search space of 403 elements:

Number of unique combinations of 2 elements: 81003
Number of unique combinations of 3 elements: 10827401
Number of unique combinations of 4 elements: 1082740100

Number of inequivalent binary ratios with max coefficient 2: 
3
Number of inequivalent binary ratios with max coefficient 4: 
11
Number of inequivalent binary ratios with max coefficient 6: 
23
Number of inequivalent binary ratios with max coefficient 8: 
43
Number of inequivalent ternary ratios with max coefficient 2: 
7
Number of inequivalent ternary ratios with max coefficient 4: 
55
Number of inequivalent ternary ratios with max coefficient 6: 
177
Number of inequivalent ternary ratios with max coefficient 8: 
433
Number of inequivalent quaternary ratios with max coefficient 2: 
15
Number of inequivalent quaternary ratios with max coefficient 4: 
239
Number of inequivalent quaternary ratios with max coefficient 6: 
1175
Number of inequivalent quaternary ratios with max coefficient 8: 
3781

Unique binary compounds (no stability constraints) with max coefficient 4: 
8.910e+05
Unique binary compounds (no stability constraints) with max coefficient 6: 
1.863e+06
Unique binary compounds (no stability constraints) with max coefficient 8: 
3.483e+06
Unique ternary compounds (no stability constraints) with max coefficient 4: 
5.955e+08
Unique ternary compounds (no stability constraints) with max coefficient 6: 
1.916e+09
Unique ternary compounds (no stability constraints) with max coefficient 8: 
4.688e+09
Unique quaternary compounds (no stability constraints) with max coefficient 4: 
2.588e+11
Unique quaternary compounds (no stability constraints) with max coefficient 6: 
1.272e+12
Unique quaternary compounds (no stability constraints) with max coefficient 8: 
4.094e+12
/Users/dan/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:18: DeprecationWarning: fractions.gcd() is deprecated. Use math.gcd() instead.