README.md

## 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](https://www.python.org/) 3.6  or later
- [NumPy](https://numpy.org/doc/stable/reference/)
- [SymPy](https://www.sympy.org/en/index.html)
- [NetworkX](https://networkx.org/)
- [pandas](https://pandas.pydata.org/)
- [itertools](https://docs.python.org/3/library/itertools.html)

### Usage
#### Make a cluster using networkx.
One can use predefined clusters in [nmf.graphs/](https://github.itap.purdue.edu/GreeleyGroup/nmf/blob/master/nmf/graphs.py)
~~~python
>>> from nmf.graphs import c_2d
>>> from nmf.util import cluster_draw
>>> cluster_draw(c_2d)
~~~
<img src="./images/Picture1.png" width="600" height="500"> </img>

#### Call the cluster class to read the graph
~~~python
>>> 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
~~~python
>>> 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
~~~
![\Large](https://latex.codecogs.com/svg.image?a_{1&space;1}^{4}&space;h_{1&space;a}^{4}&space;y_{a}^{5}&space;&plus;&space;a_{1&space;1}^{4}&space;y_{a}^{4}&space;&plus;&space;4&space;a_{1&space;1}^{3}&space;h_{1&space;a}^{3}&space;y_{a}^{4}&space;&plus;&space;4&space;a_{1&space;1}^{3}&space;y_{a}^{3}&space;&plus;&space;6&space;a_{1&space;1}^{2}&space;h_{1&space;a}^{2}&space;y_{a}^{3}&space;&plus;&space;6&space;a_{1&space;1}^{2}&space;y_{a}^{2}&space;&plus;&space;4&space;a_{1&space;1}&space;h_{1&space;a}&space;y_{a}^{2}&space;&plus;&space;4&space;a_{1&space;1}&space;y_{a}&space;&plus;&space;y_{a}&space;&plus;&space;1)
#### Get the consistency equations
~~~python
>>> 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)]
~~~
![\Large](https://latex.codecogs.com/svg.image?\inline&space;a_{1&space;1}^{4}&space;y_{a}^{4}&space;-&space;a_{1&space;1}^{3}&space;h_{1&space;a}^{3}&space;y_{a}^{4}&space;&plus;&space;3&space;a_{1&space;1}^{3}&space;y_{a}^{3}&space;-&space;3&space;a_{1&space;1}^{2}&space;h_{1&space;a}^{2}&space;y_{a}^{3}&space;&plus;&space;3&space;a_{1&space;1}^{2}&space;y_{a}^{2}&space;-&space;3&space;a_{1&space;1}&space;h_{1&space;a}&space;y_{a}^{2}&space;&plus;&space;a_{1&space;1}&space;y_{a}&space;-&space;y_{a}&space;=&space;0)

The equations generated by [cluster.consistency_equations()](https://github.itap.purdue.edu/GreeleyGroup/nmf/blob/master/nmf/nmf.py#L122) are [SymPy Equality Equations](https://docs.sympy.org/latest/modules/core.html#equality). 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
~~~python
>>> 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](https://docs.sympy.org/latest/modules/solvers/solvers.html?highlight=nsolve#sympy.solvers.solvers.nsolve) to get the coverage of species for given chemical potential and interaction energies.
~~~python
>>> 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 (θ)')
~~~
![2d cluster graph](/images/Picture2.png)