nmf python package helps to construct the partition function for a given cluster and subsequently solve the consistency equations.
One can use predefined clusters in nmf.graphs/
>>> from nmf.graphs import c_2d
>>> from nmf.util import cluster_draw
>>> cluster_draw(c_2d)
>>> from nmf import cluster
>>> clstr = cluster(c_2d,1) #the inputs are (graph , number of adsorbates)
Known Terms
Ads Symbol
a y_a
Ads Ads Bond Interaction-Symbol
a a 1 h_1_a
Unknown Terms
Ads Site Energy-Symbol
a e1 a_1_1
Calculating partition function
Calculating Consistency Equations
Eqn1: node 1 - node 0 == 0 (Site : e1, Ads : a)>>> clstr.total_expr
a_1_1**4*h_1_a**4*y_a**5 + a_1_1**4*y_a**4 + 4*a_1_1**3*h_1_a**3*y_a**4 + 4*a_1_1**3*y_a**3 + 6*a_1_1**2*h_1_a**2*y_a**3 + 6*a_1_1**2*y_a**2 + 4*a_1_1*h_1_a*y_a**2 + 4*a_1_1*y_a + y_a + 1>>> clstr.const_eqns
[Eq(a_1_1**4*y_a**4 - a_1_1**3*h_1_a**3*y_a**4 + 3*a_1_1**3*y_a**3 - 3*a_1_1**2*h_1_a**2*y_a**3 + 3*a_1_1**2*y_a**2 - 3*a_1_1*h_1_a*y_a**2 + a_1_1*y_a - y_a, 0)]The equations generated by cluster.consistency_equations() are SymPy Equality Equations. Therefore we can use SymPy Solve functionalities to find the unknowns. The unknwons are stored in list format at cluster.a_list
>>> import sympy as sp
>>> answer = sp.solve(eqns, clstr.a_list)However this takes unrealistic amount of time to solve even simple clusters. Threrfore we just use a numeric solver.
We have implemented a sympy numeric solver to get the coverage of species for given chemical potential and interaction energies.
>>> import numpy as np
>>> from matplotlib import pyplot as plt
>>> ln_mu = np.linspace(-5,45)
theta = np.zeros(ln_mu.shape)
for idx,i in enumerate(ln_mu):
params = {'y_a':np.exp(i),'h_1_a':np.exp(-8)}
theta_df = clstr.get_theta(params,maxsteps=1000)
theta[idx] = float(theta_df['e1']['a'])
plt.plot(ln_mu,theta,'-b')
plt.xlabel('Chemical Potential (Units of kT)')
plt.ylabel('Coverage (θ)')