Non-Mean-Field Cluster Calculator (nmf)

nmf python package helps to construct the partition function for a given cluster and subsequently solve the consistency equations.

Requirements

• Python 3.6 or later
• NumPy
• SymPy
• NetworkX
• pandas
• itertools

Usage

Make a cluster using networkx.

One can use predefined clusters in nmf.graphs/

>>> from nmf.graphs import c_2d
>>> from nmf.util import cluster_draw
>>> cluster_draw(c_2d)

Call the cluster class to read the graph

>>> 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)

Get the partition function

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

Get the consistency equations

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

Solve using SymPy

>>> 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.

Get theta for given chemical potential

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 (θ)')