From 5e9dc935f54622448dc41d47645fcee8fc06634f Mon Sep 17 00:00:00 2001
From: Purva Paranjape <pparanja@purdue.edu>
Date: Wed, 15 Mar 2023 10:25:49 -0400
Subject: [PATCH 1/9] Update nmf.py

---
 nmf/nmf.py | 75 ++++++++++++++++++++++++++++++++++++------------------
 1 file changed, 50 insertions(+), 25 deletions(-)

diff --git a/nmf/nmf.py b/nmf/nmf.py
index 5524c8f..7670932 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -18,11 +18,27 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
             self.names = ads_names
         else:
             self.names = [alphabet[i] for i in range(self.num_ads)]
-        self.end_sites = set([k['typ'] for i,k in self.graph.nodes(data=True)])
-        if central_site:
-            self.end_sites.remove(central_site)
-        else:
-            self.end_sites.remove("c")
+        self.site_type = set ([k['typ1'] for i,k in self.graph.nodes (data = True)])
+        print(self.site_type)
+        self.all_sites =  set ([k['typ2'] for i,k in self.graph.nodes (data = True)])
+        self.end_sites =  set ([k['typ2'] for i,k in self.graph.nodes (data = True)])
+        self.central_sites = []
+        for i,site in enumerate (self.all_sites):
+            if site.startswith ("c"):
+                self.central_sites.append(site)
+                self.end_sites.remove(site)
+        self.end_type = []
+        self.etype = []
+        for i,typ in enumerate(self.site_type):
+            e_i = []
+            for j in range(0,self.graph.number_of_nodes()):
+                if self.graph.nodes[j]['typ1'] == typ:
+                    if self.graph.nodes[j]['typ2'] in self.end_sites:
+                        e_s = self.graph.nodes[j]['typ2']
+                        e_i.append(e_s)
+                   
+            self.end_type.append(e_i)   
+        print (self.end_type)   
         self.bond_types = set([k['weight'] for i,j,k in self.graph.edges(data=True)])
         self.params = {}
         self.make_y() #Energy for adsorption
@@ -38,10 +54,17 @@ def __repr__(self):
     def make_y(self):
         #Return list of symbol y with length of list equal to number of adsorbates.
         self.y=[]
-        for i in range(0,self.num_ads):
-            y_symbol = sp.symbols(f"y_{self.names[i]}")
-            self.y.append(y_symbol)
-            self.params[f"{y_symbol}"]= y_symbol
+        self.activities=[]
+        for i in range (0,self.num_ads): 
+            y_i = {}
+            for j,site in enumerate(self.site_type):
+                y_i_symbol = sp.symbols(f"y_{site}_{i+1}", real=True)
+                y_i[f'{site}'] = y_i_symbol
+                self.activities.append(y_i_symbol)
+                self.params[f"{y_i_symbol}"]= y_i_symbol
+            self.y.append(y_i)
+            
+        print(self.y)
         return
     
     def make_h(self):
@@ -65,23 +88,23 @@ def make_h(self):
                     df[sets[1]][sets[0]]=h_dual
                     self.params[f"{h_dual}"]= h_dual
             self.h[f'{bond}']=df
+        print(self.h)
         return
-    
+     
     def make_a(self):
         '''
-        Return list of dictionaries with length of list equal to number of adsorbates.
-        Each dictionary consists of symbol for each end node type.
+        Return list of list of pandas dataframe with length of list equal to number of adsorbates.
+        Each dataframe is (# of adsorbates x # of 
         '''
-        self.a=[]
-        self.unknown_list=[]
+        self.a = {}
+        
+        columns = self.site_type
         for i in range(0,self.num_ads):
-            a_i = {}
-            for j,site in enumerate(self.end_sites):
-                a_i_symbol = sp.symbols(f"a_{i+1}_{j+1}", real=True)
-                a_i[f'{site}'] = a_i_symbol
-                self.unknown_list.append(a_i_symbol)
-                self.params[f"{a_i_symbol}"]= a_i_symbol
-            self.a.append(a_i)
+            for j,site in enumerate(self.end_type[i]):
+                df = pd.DataFrame(columns = columns)
+                for k, site in enumerate(self.site_type):
+                    a_i_sym = sp.symbols(f"a_{site}_{i+1}_{j+1}", real=True)
+                    df[self.names[i]][self.end_type[i][j]] = a_i_sym
         return
         
     def get_edge_value(self,u,v,w,G):
@@ -111,10 +134,11 @@ def make_partition_function(self):
                     for i in range(self.num_ads): #loop over number of adsorbates
                         for j in lists[i]: #loop over sites
                             G.nodes[j]["occ"]=self.names[i]
-                            th *= self.y[i]
-                            for k in self.end_sites: #loop over end types
-                                if G.nodes[j]["typ"]==k:
-                                    th *= self.a[i][k]
+                            for k in self.site_type:
+                                    th *= self.y[i][k]
+                                    for k, l in self.end_typ: #loop over end types
+                                        if G.nodes[j]["typ"].endswith(l):
+                                            th *= self.a[i][l]
                     for u,v,w in G.edges(data=True):
                         th *= self.get_edge_value(u,v,w,G)
                     self.total_expr += th
@@ -229,3 +253,4 @@ def pprint(self):
     
     def draw(self):
         return cluster_draw(self.graph)
+

From 47c4ed6c0a5969fe8a3c36e95fcf7b5c0d8eec43 Mon Sep 17 00:00:00 2001
From: Purva Paranjape <pparanja@purdue.edu>
Date: Wed, 15 Mar 2023 17:05:37 -0400
Subject: [PATCH 2/9] Update nmf.py

updated activities for different binding modes
---
 nmf/nmf.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/nmf/nmf.py b/nmf/nmf.py
index 7670932..6a0bad4 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -98,13 +98,14 @@ def make_a(self):
         '''
         self.a = {}
         
-        columns = self.site_type
-        for i in range(0,self.num_ads):
-            for j,site in enumerate(self.end_type[i]):
-                df = pd.DataFrame(columns = columns)
-                for k, site in enumerate(self.site_type):
-                    a_i_sym = sp.symbols(f"a_{site}_{i+1}_{j+1}", real=True)
+        columns = self.num_ads
+        for i,site in enumerate(self.site_type):
+            for j in range(self.num_ads):
+                df = pd.DataFrame(index = self.end_type[i],columns = columns)
+                for k, site in enumerate(self.end_type[i]):
+                    a_i_sym = sp.symbols(f"a_{site}_{j+1}_{k+1}", real=True)
                     df[self.names[i]][self.end_type[i][j]] = a_i_sym
+        a.append(a_i_sym)
         return
         
     def get_edge_value(self,u,v,w,G):
@@ -137,7 +138,7 @@ def make_partition_function(self):
                             for k in self.site_type:
                                     th *= self.y[i][k]
                                     for k, l in self.end_typ: #loop over end types
-                                        if G.nodes[j]["typ"].endswith(l):
+                                        if G.nodes[j][k]:
                                             th *= self.a[i][l]
                     for u,v,w in G.edges(data=True):
                         th *= self.get_edge_value(u,v,w,G)

From 68b54250ea48c83fde0f8c86fea13d39798583b9 Mon Sep 17 00:00:00 2001
From: Purva Paranjape <pparanja@purdue.edu>
Date: Wed, 15 Mar 2023 21:34:16 -0400
Subject: [PATCH 3/9] Update nmf.py

---
 nmf/nmf.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/nmf/nmf.py b/nmf/nmf.py
index 6a0bad4..e9f8e6f 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -98,14 +98,16 @@ def make_a(self):
         '''
         self.a = {}
         
-        columns = self.num_ads
+        columns = [self.names[i] for i in range(self.num_ads)]
         for i,site in enumerate(self.site_type):
-            for j in range(self.num_ads):
-                df = pd.DataFrame(index = self.end_type[i],columns = columns)
-                for k, site in enumerate(self.end_type[i]):
+            df = pd.DataFrame(index = self.end_type[i],columns = columns)
+            print(df)
+            for k, end in enumerate(self.end_type[i]):
+                for j,ads in enumerate(columns):
+                    print(site,end,ads)
                     a_i_sym = sp.symbols(f"a_{site}_{j+1}_{k+1}", real=True)
-                    df[self.names[i]][self.end_type[i][j]] = a_i_sym
-        a.append(a_i_sym)
+                    df[ads][end] = a_i_sym
+            self.a["site"] = df
         return
         
     def get_edge_value(self,u,v,w,G):

From fecf763f0cdc0aba63c15c64d5d493a52dc829ff Mon Sep 17 00:00:00 2001
From: Purva Paranjape <pparanja@purdue.edu>
Date: Thu, 16 Mar 2023 16:54:37 -0400
Subject: [PATCH 4/9] Update nmf.py

---
 nmf/nmf.py | 39 ++++++++++++++++++++++-----------------
 1 file changed, 22 insertions(+), 17 deletions(-)

diff --git a/nmf/nmf.py b/nmf/nmf.py
index e9f8e6f..ebb1c93 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -46,7 +46,7 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
         self.make_a() #Energy of end atoms
         self.pprint() #Print details of class
         self.make_partition_function() 
-        self.make_consistency_equations()
+       # self.make_consistency_equations()
         
     def __repr__(self):
         return 'Cluster'
@@ -97,7 +97,7 @@ def make_a(self):
         Each dataframe is (# of adsorbates x # of 
         '''
         self.a = {}
-        
+        self.unknown_list = []
         columns = [self.names[i] for i in range(self.num_ads)]
         for i,site in enumerate(self.site_type):
             df = pd.DataFrame(index = self.end_type[i],columns = columns)
@@ -105,9 +105,11 @@ def make_a(self):
             for k, end in enumerate(self.end_type[i]):
                 for j,ads in enumerate(columns):
                     print(site,end,ads)
-                    a_i_sym = sp.symbols(f"a_{site}_{j+1}_{k+1}", real=True)
+                    a_i_sym = sp.symbols(f"a_{site}_{ads}_{k+1}", real=True)
+                    self.unknown_list.append(a_i_sym)
                     df[ads][end] = a_i_sym
-            self.a["site"] = df
+            self.a[f'{site}'] = df
+        print(self.a)
         return
         
     def get_edge_value(self,u,v,w,G):
@@ -137,18 +139,19 @@ def make_partition_function(self):
                     for i in range(self.num_ads): #loop over number of adsorbates
                         for j in lists[i]: #loop over sites
                             G.nodes[j]["occ"]=self.names[i]
-                            for k in self.site_type:
-                                    th *= self.y[i][k]
-                                    for k, l in self.end_typ: #loop over end types
-                                        if G.nodes[j][k]:
-                                            th *= self.a[i][l]
+                            for k, site in enumerate(self.site_type):
+                                if G.nodes[j]['typ1'] == site:
+                                    th *= self.y[i][site]
+                                    for l,end in enumerate(self.end_type[k]): #loop over end types
+                                        if G.nodes[j]['typ2'] == end:
+                                            th *= self.a[site][self.names[i]][end]
                     for u,v,w in G.edges(data=True):
                         th *= self.get_edge_value(u,v,w,G)
                     self.total_expr += th
                     self.complete_site_list.append({'graph':G.copy(),'expr':th,})
                     for i in G.nodes():
                         G.nodes[i]["occ"]='-'
-        #print(f"z={self.total_expr}")
+        print(f"z={self.total_expr}")
         return self.total_expr
     
     def make_consistency_equations(self):
@@ -177,7 +180,7 @@ def make_consistency_equations(self):
         for ads in range(self.num_ads):
             for site in self.end_sites:
                 for node,data in self.graph.nodes(data=True):
-                    if data['typ'] == site:
+                    if data['typ2'] == site:
                         print(f"Eqn{idx}: node {node} - node 0 == 0 (Site : {site}, Ads : {self.names[ads]})")
                         self.const_eqns.append(sp.Eq(self.node_eqns[node][self.names[ads]]-self.node_eqns[0][self.names[ads]],0))
                         self.const_eqns_expr.append(self.node_eqns[node][self.names[ads]]-self.node_eqns[0][self.names[ads]])
@@ -239,7 +242,7 @@ def pprint(self):
         print("Ads\tSymbol")
         for idx,i in enumerate(self.y):
             print(f"{self.names[idx]}\t{i}")
-        print("\nAds\tAds\tBond\tInteraction-Symbol")
+        print("\nAdsorbate\tAdsorbate\tBond\tInteraction Energy")
         for bond in self.h:
             if bond!='0':
                 for ads1 in self.h[bond]:
@@ -248,12 +251,14 @@ def pprint(self):
                         if h_symbol in self.unknown_list:
                             print(f"{h_symbol} is unknown")
         print("\nUnknown Terms")
-        print("Ads\tSite\tEnergy-Symbol")
-        for idx,ads1 in enumerate(self.a):
-            for key,value in ads1.items():
-                print(f"{self.names[idx]}\t{key}\t{value}")
+        print("Site\tAdsorbate\tEdge\tCorrection Field")
+        for site in self.a:
+            for ads in self.a[site]:
+                for end, a_sym in self.a[site][ads].iteritems():
+                    print(f"{site}\t{ads}\t{end}\t{a_sym}")
+        
         return
-    
+                                        
     def draw(self):
         return cluster_draw(self.graph)
 

From 3fab8c19bbc63b5502c6365c97cff00dfc29d72e Mon Sep 17 00:00:00 2001
From: Purva Paranjape <pparanja@purdue.edu>
Date: Fri, 7 Apr 2023 15:48:54 -0400
Subject: [PATCH 5/9] Update graphs.py

modified graphs to include typ1 & typ2
---
 nmf/graphs.py | 48 ++++++++++++++++++++++++------------------------
 1 file changed, 24 insertions(+), 24 deletions(-)

diff --git a/nmf/graphs.py b/nmf/graphs.py
index 6a3034b..e520e73 100644
--- a/nmf/graphs.py
+++ b/nmf/graphs.py
@@ -8,18 +8,18 @@
 
 #1d Graph square symmetry
 c_1d = nx.Graph()
-c_1d.add_node(0,typ="c")
-c_1d.add_node(1,typ="e1")
-c_1d.add_node(2,typ="e1")
+c_1d.add_node(0,typ1 = "t", typ2 ="c1")
+c_1d.add_node(1,typ1 = "t", typ2 ="e1")
+c_1d.add_node(2,typ1 = "t", typ2 ="e2")
 c_1d.add_weighted_edges_from([(0, 2, 1), (0, 1, 1)])
 
 #2d Graph square symmetry
 c_2d = nx.Graph()
-c_2d.add_node(0,typ="c")
-c_2d.add_node(1,typ="e1")
-c_2d.add_node(2,typ="e1")
-c_2d.add_node(3,typ="e1")
-c_2d.add_node(4,typ="e1")
+c_2d.add_node(0,typ1 = "t", typ2 ="c1")
+c_2d.add_node(1,typ1 = "t", typ2 ="e1")
+c_2d.add_node(2,typ1 = "t", typ2 ="e1")
+c_2d.add_node(3,typ1 = "t", typ2 ="e1")
+c_2d.add_node(4,typ1 = "t", typ2 ="e1")
 c_2d.add_weighted_edges_from([(0, 2, 1), (0, 1, 1), (0, 3, 1), (0, 4, 1)])
 
 #2d Graph with edge nodes connected to each other
@@ -28,13 +28,13 @@
 
 #2d Graph Hexagonal symmetry
 c_2d_hex = nx.Graph()
-c_2d_hex.add_node(0,typ="c")
-c_2d_hex.add_node(1,typ="e1")
-c_2d_hex.add_node(2,typ="e1")
-c_2d_hex.add_node(3,typ="e1")
-c_2d_hex.add_node(4,typ="e1")
-c_2d_hex.add_node(5,typ="e1")
-c_2d_hex.add_node(6,typ="e1")
+c_2d_hex.add_node(0,typ1 = "t", typ2 ="c1")
+c_2d_hex.add_node(1,typ1 = "t", typ2 ="e1")
+c_2d_hex.add_node(2,typ1 = "t", typ2 ="e1")
+c_2d_hex.add_node(3,typ1 = "t", typ2 ="e1")
+c_2d_hex.add_node(4,typ1 = "t", typ2 ="e1")
+c_2d_hex.add_node(5,typ1 = "t", typ2 ="e1")
+c_2d_hex.add_node(6,typ1 = "t", typ2 ="e1")
 c_2d_hex.add_weighted_edges_from([(0, 2, 1), (0, 1, 1), (0, 3, 1), (0, 4, 1), (0, 5, 1), (0, 6, 1)])
 
 #2d Hexagonal Graph with edge nodes connected to each other
@@ -44,15 +44,15 @@
 #2d Graph with asymmetric bridge sites
 custom_1 = nx.Graph()
 custom_1 = nx.Graph()
-custom_1.add_node(0,typ="c")
-custom_1.add_node(2,typ="e1")
-custom_1.add_node(4,typ="e1")
-custom_1.add_node(6,typ="e1")
-custom_1.add_node(8,typ="e1")
-custom_1.add_node(1,typ="e2")
-custom_1.add_node(5,typ="e2")
-custom_1.add_node(3,typ="e3")
-custom_1.add_node(7,typ="e3")
+custom_1.add_node(0,typ1 = "b",typ2="c1")
+custom_1.add_node(2,typ1 = "b",typ2="e1")
+custom_1.add_node(4,typ1 = "b",typ2="e1")
+custom_1.add_node(6,typ1 = "b",typ2="e1")
+custom_1.add_node(8,typ1 = "b",typ2="e1")
+custom_1.add_node(1,typ1 = "b",typ2="e2")
+custom_1.add_node(5,typ1 = "b",typ2="e2")
+custom_1.add_node(3,typ1 = "b",typ2="e3")
+custom_1.add_node(7,typ1 = "b",typ2="e3")
 custom_1.add_weighted_edges_from([ (0, 1, 2), (0, 7, 1), (0, 3, 1), (0, 5, 2), 
                             (8, 2, 1), (2, 4, 2), (6, 4, 1), (6, 8, 2),
                             (2, 1, 0), (2, 3, 0), (3, 4, 0), (4, 5, 0),

From 34323528e7602d9dddaffa0c017f729dfa9cfdb6 Mon Sep 17 00:00:00 2001
From: kkjsawant <sawantk@purdue.edu>
Date: Wed, 19 Apr 2023 16:10:59 -0400
Subject: [PATCH 6/9] Fixed minor issues

---
 nmf/graphs.py |   2 +-
 nmf/nmf.py    | 107 +++++++++++++++++++++++---------------------------
 2 files changed, 51 insertions(+), 58 deletions(-)

diff --git a/nmf/graphs.py b/nmf/graphs.py
index e520e73..ec56a7e 100644
--- a/nmf/graphs.py
+++ b/nmf/graphs.py
@@ -10,7 +10,7 @@
 c_1d = nx.Graph()
 c_1d.add_node(0,typ1 = "t", typ2 ="c1")
 c_1d.add_node(1,typ1 = "t", typ2 ="e1")
-c_1d.add_node(2,typ1 = "t", typ2 ="e2")
+c_1d.add_node(2,typ1 = "t", typ2 ="e1")
 c_1d.add_weighted_edges_from([(0, 2, 1), (0, 1, 1)])
 
 #2d Graph square symmetry
diff --git a/nmf/nmf.py b/nmf/nmf.py
index ebb1c93..7d2fdc8 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -14,21 +14,24 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
         for node in self.graph:
             self.graph.nodes[node]['occ'] = '-'
         self.num_ads = num_ads
+        
+        #If adsorbates names are given they are stored as such, otherwise we use alphabets starting from a,b,...
         if ads_names:
             self.names = ads_names
         else:
             self.names = [alphabet[i] for i in range(self.num_ads)]
-        self.site_type = set ([k['typ1'] for i,k in self.graph.nodes (data = True)])
-        print(self.site_type)
-        self.all_sites =  set ([k['typ2'] for i,k in self.graph.nodes (data = True)])
-        self.end_sites =  set ([k['typ2'] for i,k in self.graph.nodes (data = True)])
+            
+        #The sites are seperated as central or end depending on how they are assigned.    
+        self.site_type = set([k['typ1'] for i,k in self.graph.nodes (data = True)])
+        self.all_sites = set([k['typ2'] for i,k in self.graph.nodes (data = True)])
+        self.end_sites = set([k['typ2'] for i,k in self.graph.nodes (data = True)])
         self.central_sites = []
         for i,site in enumerate (self.all_sites):
-            if site.startswith ("c"):
+            if site.startswith ("c"): #If the typ2 attribute starts with c it is stored as a central site
                 self.central_sites.append(site)
                 self.end_sites.remove(site)
-        self.end_type = []
-        self.etype = []
+                
+        self.end_type = {}
         for i,typ in enumerate(self.site_type):
             e_i = []
             for j in range(0,self.graph.number_of_nodes()):
@@ -37,8 +40,7 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
                         e_s = self.graph.nodes[j]['typ2']
                         e_i.append(e_s)
                    
-            self.end_type.append(e_i)   
-        print (self.end_type)   
+            self.end_type[f"{typ}"]=set(e_i)    
         self.bond_types = set([k['weight'] for i,j,k in self.graph.edges(data=True)])
         self.params = {}
         self.make_y() #Energy for adsorption
@@ -46,25 +48,23 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
         self.make_a() #Energy of end atoms
         self.pprint() #Print details of class
         self.make_partition_function() 
-       # self.make_consistency_equations()
+        self.make_consistency_equations()
         
     def __repr__(self):
         return 'Cluster'
     
     def make_y(self):
-        #Return list of symbol y with length of list equal to number of adsorbates.
-        self.y=[]
+        '''
+        Returns a list of lists. Where first index is name of adsorbate and the sec
+        '''
+        self.y=pd.DataFrame(index=self.names,columns=self.site_type)
         self.activities=[]
-        for i in range (0,self.num_ads): 
-            y_i = {}
-            for j,site in enumerate(self.site_type):
-                y_i_symbol = sp.symbols(f"y_{site}_{i+1}", real=True)
-                y_i[f'{site}'] = y_i_symbol
+        for ads in self.names:
+            for site in self.site_type:
+                y_i_symbol = sp.symbols(f"y_{site}_{ads}", real=True)
+                self.y[f'{site}'][f'{ads}'] = y_i_symbol
                 self.activities.append(y_i_symbol)
                 self.params[f"{y_i_symbol}"]= y_i_symbol
-            self.y.append(y_i)
-            
-        print(self.y)
         return
     
     def make_h(self):
@@ -88,28 +88,24 @@ def make_h(self):
                     df[sets[1]][sets[0]]=h_dual
                     self.params[f"{h_dual}"]= h_dual
             self.h[f'{bond}']=df
-        print(self.h)
         return
      
     def make_a(self):
         '''
-        Return list of list of pandas dataframe with length of list equal to number of adsorbates.
-        Each dataframe is (# of adsorbates x # of 
+        Return dictionary of pandas dataframe.
+        Each dataframe is (# of adsorbates , # of end sites)
+        Each key in ditionary is site type
         '''
         self.a = {}
         self.unknown_list = []
-        columns = [self.names[i] for i in range(self.num_ads)]
-        for i,site in enumerate(self.site_type):
-            df = pd.DataFrame(index = self.end_type[i],columns = columns)
-            print(df)
-            for k, end in enumerate(self.end_type[i]):
-                for j,ads in enumerate(columns):
-                    print(site,end,ads)
-                    a_i_sym = sp.symbols(f"a_{site}_{ads}_{k+1}", real=True)
+        for etype in self.end_type.keys():
+            df = pd.DataFrame(index = self.end_type[f"{etype}"],columns = self.names)
+            for end in self.end_type[f"{etype}"]:
+                for ads in self.names:
+                    a_i_sym = sp.symbols(f"a_{etype}_{ads}_{end}", real=True)
                     self.unknown_list.append(a_i_sym)
                     df[ads][end] = a_i_sym
-            self.a[f'{site}'] = df
-        print(self.a)
+            self.a[f'{etype}'] = df
         return
         
     def get_edge_value(self,u,v,w,G):
@@ -136,15 +132,15 @@ def make_partition_function(self):
                 all_lists = divide_adsorbate(self.num_ads,sets) #Divide sites among adsorbates
                 for lists in all_lists: #loop over all the possibilities
                     th=1 #Intitalize site term
-                    for i in range(self.num_ads): #loop over number of adsorbates
+                    for i,ads in enumerate(self.names): #loop over number of adsorbates
                         for j in lists[i]: #loop over sites
-                            G.nodes[j]["occ"]=self.names[i]
-                            for k, site in enumerate(self.site_type):
+                            G.nodes[j]["occ"]=ads
+                            for site in self.site_type:
                                 if G.nodes[j]['typ1'] == site:
-                                    th *= self.y[i][site]
-                                    for l,end in enumerate(self.end_type[k]): #loop over end types
+                                    th *= self.y[site][ads]
+                                    for end in self.end_sites: #loop over end types
                                         if G.nodes[j]['typ2'] == end:
-                                            th *= self.a[site][self.names[i]][end]
+                                            th *= self.a[site][ads][end]
                     for u,v,w in G.edges(data=True):
                         th *= self.get_edge_value(u,v,w,G)
                     self.total_expr += th
@@ -193,13 +189,8 @@ def make_consistency_equations(self):
         self.ln = 0
         return
     
-    def get_theta(self,params, initial_guess=None, maxsteps=50, check_params=True, verify=False):
-        if check_params:
-            key_list = list(params.keys())
-            for i in self.y:
-                assert str(i) in key_list, "Missing Input Values"
+    def get_theta(self, params, initial_guess=None, maxsteps=50, check_params=True, verify=False):
 
-        #sub_eqns = [None] * len(self.const_eqns)
         sub_eqns_expr = [None] * len(self.const_eqns)
         initial_guess1 = [0.5] * len(self.unknown_list)
         if initial_guess:
@@ -239,24 +230,26 @@ def add_unknown(self,symbol,eqn):
     def pprint(self):
         #Print all the symbols used.
         print("Known Terms")
-        print("Ads\tSymbol")
-        for idx,i in enumerate(self.y):
-            print(f"{self.names[idx]}\t{i}")
-        print("\nAdsorbate\tAdsorbate\tBond\tInteraction Energy")
+        print("Site\tAds\tSymbol")
+        for site in self.site_type:
+            for ads in self.names:
+                y = self.y[f"{site}"][f"{ads}"]
+                print(f"{site}\t{ads}\t{y}")
+        print("\nAds\tAds\tBond\tInteraction Energy")
         for bond in self.h:
             if bond!='0':
                 for ads1 in self.h[bond]:
                     for ads2, h_symbol in self.h[bond][ads1].iteritems():
                         print(f"{ads1}\t{ads2}\t{bond}\t{h_symbol}")
-                        if h_symbol in self.unknown_list:
-                            print(f"{h_symbol} is unknown")
-        print("\nUnknown Terms")
-        print("Site\tAdsorbate\tEdge\tCorrection Field")
-        for site in self.a:
-            for ads in self.a[site]:
-                for end, a_sym in self.a[site][ads].iteritems():
-                    print(f"{site}\t{ads}\t{end}\t{a_sym}")
         
+        print("Unknown Terms")
+        print("Site\tAds\tEndSite\tSymbol")
+        for site in self.end_type.keys():
+            for end in self.end_type[f"{site}"]:
+                for ads in self.names:
+                    a = self.a[f"{site}"][f"{ads}"][f"{end}"]
+                    print(f"{site}\t{ads}\t{end}\t{a}")
+       
         return
                                         
     def draw(self):

From 24f234695f086822370576713e74f7b4fa94fee1 Mon Sep 17 00:00:00 2001
From: kkjsawant <sawantk@purdue.edu>
Date: Tue, 2 May 2023 14:38:52 -0400
Subject: [PATCH 7/9] added multidentate

---
 nmf/nmf.py   | 96 +++++++++++++++++++++++++++++++++++++++++-----------
 nmf/untitled |  0
 nmf/util.py  | 11 ++++++
 3 files changed, 88 insertions(+), 19 deletions(-)
 create mode 100644 nmf/untitled

diff --git a/nmf/nmf.py b/nmf/nmf.py
index 7d2fdc8..9420557 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -3,11 +3,11 @@
 import sympy as sp
 import pandas as pd
 from itertools import combinations
-from nmf.util import alphabet,divide_adsorbate,cluster_draw
+from nmf.util import alphabet,divide_adsorbate,cluster_draw,divide_stoich
 from scipy.optimize import fsolve
 
 class cluster:
-    def __init__(self,graph,num_ads,ads_names=None,central_site=None):
+    def __init__(self,graph,num_ads,ads_names=None,ads_stoich=None,central_site=None):
         
         assert type(graph).__name__ == 'Graph',"Not a Graph"
         self.graph = graph
@@ -20,6 +20,13 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
             self.names = ads_names
         else:
             self.names = [alphabet[i] for i in range(self.num_ads)]
+        
+        #To deal with dual sites
+        if ads_stoich:
+            assert len(self.names) == len(ads_stoich), "Length of ads_stoich doesnt match with number of adsorbates"
+            self.ads_stoich = ads_stoich
+        else:
+            self.ads_stoich = [1] * int(self.num_ads)
             
         #The sites are seperated as central or end depending on how they are assigned.    
         self.site_type = set([k['typ1'] for i,k in self.graph.nodes (data = True)])
@@ -48,7 +55,7 @@ def __init__(self,graph,num_ads,ads_names=None,central_site=None):
         self.make_a() #Energy of end atoms
         self.pprint() #Print details of class
         self.make_partition_function() 
-        self.make_consistency_equations()
+        #self.make_consistency_equations()
         
     def __repr__(self):
         return 'Cluster'
@@ -88,6 +95,11 @@ def make_h(self):
                     df[sets[1]][sets[0]]=h_dual
                     self.params[f"{h_dual}"]= h_dual
             self.h[f'{bond}']=df
+            
+        df = pd.DataFrame(index=index,columns=index)
+        for i in range(self.num_ads):
+            df[self.names[i]][self.names[i]]=1
+        self.h['-1']=df
         return
      
     def make_a(self):
@@ -127,26 +139,71 @@ def make_partition_function(self):
         self.complete_site_list=[]
         G = self.graph.copy()
         self.total_expr = 1 #Intitialize the expression (1:bare site)
-        for occ in range(1,len(G)+1): #Loop over number of sites to be filled
-            for sets in combinations(G,occ): #loop over combinations for given number of site filling
+        for occ in range(1,len(self.graph)+1): #Loop over number of sites to be filled
+            for sets in combinations(self.graph,occ): #loop over combinations for given number of site filling
                 all_lists = divide_adsorbate(self.num_ads,sets) #Divide sites among adsorbates
                 for lists in all_lists: #loop over all the possibilities
+                    G = self.graph.copy()
+                    new_graphs = []
                     th=1 #Intitalize site term
                     for i,ads in enumerate(self.names): #loop over number of adsorbates
-                        for j in lists[i]: #loop over sites
-                            G.nodes[j]["occ"]=ads
-                            for site in self.site_type:
-                                if G.nodes[j]['typ1'] == site:
-                                    th *= self.y[site][ads]
-                                    for end in self.end_sites: #loop over end types
-                                        if G.nodes[j]['typ2'] == end:
-                                            th *= self.a[site][ads][end]
-                    for u,v,w in G.edges(data=True):
-                        th *= self.get_edge_value(u,v,w,G)
+                        if self.ads_stoich[i] == 1:
+                            for j in lists[i]: #loop over sites
+                                G.nodes[j]["occ"]=ads
+                                site = G.nodes[j]['typ1']
+                                th *= self.y[site][ads]
+                                end = G.nodes[j]['typ2']
+                                if end in self.end_sites:
+                                    th *= self.a[site][ads][end]
+                                                
+                        else:
+                            if len(lists[i])%self.ads_stoich[i] ==0:
+                                new_expr = 0
+                                num_pairs = int(len(lists[i])/self.ads_stoich[i])
+                                new_list = divide_stoich(num_pairs,lists[i],stoich=self.ads_stoich[i])
+                                for poss_solution in new_list:
+                                    flag = False
+                                    n_G = G.copy()
+                                    new_th=1
+                                    for set_pair in poss_solution:
+                                        for pair in combinations(set_pair,2):
+                                            if n_G.nodes[pair[0]]["occ"]=='-' and n_G.nodes[pair[1]]["occ"]=='-' and n_G.has_edge(pair[0],pair[1]):
+                                                n_G.nodes[pair[0]]["occ"]=ads
+                                                n_G.nodes[pair[1]]["occ"]=ads
+                                                n_G[pair[0]][pair[1]]['weight']=-1
+                                                site1, site2 = n_G.nodes[pair[0]]['typ1'], n_G.nodes[pair[1]]['typ1']
+                                                new_th *= (self.y[site1][ads]/2 + self.y[site2][ads]/2)
+                                                end1, end2 = n_G.nodes[pair[0]]['typ2'], n_G.nodes[pair[1]]['typ2']
+                                                if end1 in self.end_sites:
+                                                    new_th *= self.a[site1][ads][end1]
+                                                if end2 in self.end_sites:
+                                                    new_th *= self.a[site2][ads][end2]
+                                                            
+                                            else:
+                                                new_th*=0
+                                                flag=True
+                                                break
+                                    if flag:
+                                        continue
+                                    else:
+                                        new_expr += new_th
+                                        valid_G = n_G.copy()
+                                
+                                if new_expr:
+                                    G = valid_G
+                                    th*= new_expr   
+                                else:
+                                    th=0
+                            else:
+                                th=0
+                                
+                    if th!=0 and th!=1:
+                        for u,v,w in G.edges(data=True):
+                            th *= self.get_edge_value(u,v,w,G)
+                    else:
+                        th=0
                     self.total_expr += th
                     self.complete_site_list.append({'graph':G.copy(),'expr':th,})
-                    for i in G.nodes():
-                        G.nodes[i]["occ"]='-'
         print(f"z={self.total_expr}")
         return self.total_expr
     
@@ -235,14 +292,15 @@ def pprint(self):
             for ads in self.names:
                 y = self.y[f"{site}"][f"{ads}"]
                 print(f"{site}\t{ads}\t{y}")
+                
         print("\nAds\tAds\tBond\tInteraction Energy")
         for bond in self.h:
-            if bond!='0':
+            if bond not in ['0','-1']:
                 for ads1 in self.h[bond]:
                     for ads2, h_symbol in self.h[bond][ads1].iteritems():
                         print(f"{ads1}\t{ads2}\t{bond}\t{h_symbol}")
         
-        print("Unknown Terms")
+        print("\nUnknown Terms")
         print("Site\tAds\tEndSite\tSymbol")
         for site in self.end_type.keys():
             for end in self.end_type[f"{site}"]:
diff --git a/nmf/untitled b/nmf/untitled
new file mode 100644
index 0000000..e69de29
diff --git a/nmf/util.py b/nmf/util.py
index d26b3c5..05c14a9 100644
--- a/nmf/util.py
+++ b/nmf/util.py
@@ -89,6 +89,17 @@ def divide_adsorbate(num_ads,positions):
         S.append(combination)
     return S
 
+def divide_stoich(num_pairs,lists,stoich=2):
+    solution = []
+    list2 = [list(set1) for set1 in combinations(lists,stoich)]
+    for set2 in combinations(list2,num_pairs):
+        sum1 = []
+        for l in set2:
+            sum1 += l
+        if len(sum1)==len(set(sum1)):
+            solution.append(set2)
+    return solution
+    
 def cluster_draw(cluster):
     labels = dict()
     pos = nx.spring_layout(cluster,seed=3)

From c558559961f8c67335a13f6370e046eea7825cfa Mon Sep 17 00:00:00 2001
From: kkjsawant <sawantk@purdue.edu>
Date: Wed, 3 May 2023 15:11:33 -0400
Subject: [PATCH 8/9] fixed multidentate

---
 nmf/nmf.py  | 83 ++++++++++++++++++++++++++++-------------------------
 nmf/util.py |  2 +-
 2 files changed, 45 insertions(+), 40 deletions(-)

diff --git a/nmf/nmf.py b/nmf/nmf.py
index 9420557..0d907a3 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -158,46 +158,49 @@ def make_partition_function(self):
                                                 
                         else:
                             if len(lists[i])%self.ads_stoich[i] ==0:
-                                new_expr = 0
-                                num_pairs = int(len(lists[i])/self.ads_stoich[i])
-                                new_list = divide_stoich(num_pairs,lists[i],stoich=self.ads_stoich[i])
-                                for poss_solution in new_list:
-                                    flag = False
-                                    n_G = G.copy()
-                                    new_th=1
-                                    for set_pair in poss_solution:
-                                        for pair in combinations(set_pair,2):
-                                            if n_G.nodes[pair[0]]["occ"]=='-' and n_G.nodes[pair[1]]["occ"]=='-' and n_G.has_edge(pair[0],pair[1]):
-                                                n_G.nodes[pair[0]]["occ"]=ads
-                                                n_G.nodes[pair[1]]["occ"]=ads
-                                                n_G[pair[0]][pair[1]]['weight']=-1
-                                                site1, site2 = n_G.nodes[pair[0]]['typ1'], n_G.nodes[pair[1]]['typ1']
-                                                new_th *= (self.y[site1][ads]/2 + self.y[site2][ads]/2)
-                                                end1, end2 = n_G.nodes[pair[0]]['typ2'], n_G.nodes[pair[1]]['typ2']
-                                                if end1 in self.end_sites:
-                                                    new_th *= self.a[site1][ads][end1]
-                                                if end2 in self.end_sites:
-                                                    new_th *= self.a[site2][ads][end2]
-                                                            
-                                            else:
-                                                new_th*=0
-                                                flag=True
-                                                break
-                                    if flag:
-                                        continue
+                                if lists[i]:
+                                    new_expr = 0
+                                    num_pairs = int(len(lists[i])/self.ads_stoich[i])
+                                    new_list = divide_stoich(num_pairs,lists[i],stoich=self.ads_stoich[i])
+                                    for poss_solution in new_list:
+                                        flag = False
+                                        n_G = G.copy()
+                                        new_th=1
+                                        for set_pair in poss_solution:
+                                            for pair in combinations(set_pair,2):
+                                                if n_G.has_edge(pair[0],pair[1]):
+                                                    n_G.nodes[pair[0]]["occ"]=ads
+                                                    n_G.nodes[pair[1]]["occ"]=ads
+                                                    n_G[pair[0]][pair[1]]['weight']=-1                
+                                                else:
+                                                    new_th*=0
+                                                    flag=True
+                                                    break
+
+                                            if not flag:
+                                                tmp = 0
+                                                for elem in set_pair:
+                                                    site1 = n_G.nodes[elem]['typ1']
+                                                    tmp += self.y[site1][ads]/self.ads_stoich[i]
+                                                    end1 = n_G.nodes[elem]['typ2']
+                                                    if end1 in self.end_sites:
+                                                        new_th *= self.a[site1][ads][end1]
+                                                new_th *= tmp
+                                        if flag:
+                                            continue
+                                        else:
+                                            new_expr += new_th
+                                            valid_G = n_G.copy()
+
+                                    if new_expr:
+                                        G = valid_G
+                                        th*= new_expr   
                                     else:
-                                        new_expr += new_th
-                                        valid_G = n_G.copy()
-                                
-                                if new_expr:
-                                    G = valid_G
-                                    th*= new_expr   
-                                else:
-                                    th=0
+                                        th=0
                             else:
                                 th=0
                                 
-                    if th!=0 and th!=1:
+                    if th not in [0,1]:
                         for u,v,w in G.edges(data=True):
                             th *= self.get_edge_value(u,v,w,G)
                     else:
@@ -286,12 +289,14 @@ def add_unknown(self,symbol,eqn):
     
     def pprint(self):
         #Print all the symbols used.
-        print("Known Terms")
-        print("Site\tAds\tSymbol")
+        print(f"The cluster has {len(self.graph)} adsorption sites and {self.num_ads} adsorbates")
+
+        print("\nKnown Terms")
+        print("Site\tAds\tSymbol\tStoich")
         for site in self.site_type:
             for ads in self.names:
                 y = self.y[f"{site}"][f"{ads}"]
-                print(f"{site}\t{ads}\t{y}")
+                print(f"{site}\t{ads}\t{y}\t{self.ads_stoich[self.names.index(ads)]}")
                 
         print("\nAds\tAds\tBond\tInteraction Energy")
         for bond in self.h:
diff --git a/nmf/util.py b/nmf/util.py
index 05c14a9..0de37b5 100644
--- a/nmf/util.py
+++ b/nmf/util.py
@@ -104,7 +104,7 @@ def cluster_draw(cluster):
     labels = dict()
     pos = nx.spring_layout(cluster,seed=3)
     for node,data in cluster.nodes(data=True):
-        labels[node] = f"{node},{data['typ']}"
+        labels[node] = f"{node},{data['typ1']},{data['typ2']}"
     options = {"edgecolors": "tab:gray", "node_size": 800, "alpha": 0.9}
     nx.draw_networkx_nodes(cluster, pos, node_color="tab:red",**options)
     nx.draw_networkx_edges(cluster, pos, width=1.0, alpha=0.5, edge_color="tab:red")

From 71ab8785b6ab9112033ec2712da7056e1c4ce000 Mon Sep 17 00:00:00 2001
From: "Paranjape, Purva Harshad" <pparanja@purdue.edu>
Date: Mon, 8 May 2023 17:07:26 -0400
Subject: [PATCH 9/9] get_theta & graphs

---
 nmf/graphs.py | 47 +++++++++++++++++++++++++++++++++++++++++++++++
 nmf/nmf.py    | 12 +++++++++---
 2 files changed, 56 insertions(+), 3 deletions(-)

diff --git a/nmf/graphs.py b/nmf/graphs.py
index ec56a7e..c6ca0e3 100644
--- a/nmf/graphs.py
+++ b/nmf/graphs.py
@@ -58,3 +58,50 @@
                             (2, 1, 0), (2, 3, 0), (3, 4, 0), (4, 5, 0),
                             (5, 6, 0), (6, 7, 0), (7, 8, 0), (8, 1, 0), 
                             (2, 0, 0), (4, 0, 0), (6, 0, 0), (8, 0, 0)])
+
+#2d graph for bidentates
+c_2d_9 = nx.Graph()
+c_2d_9.add_node(0,typ1 = "t", typ2 ="c1")
+c_2d_9.add_node(1,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(2,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(3,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(4,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(5,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(6,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(7,typ1 = "t", typ2 ="e1")
+c_2d_9.add_node(8,typ1 = "t", typ2 ="e1")
+c_2d_9.add_weighted_edges_from([(0,1,1),(0,3,1),(0,5,1),(0,7,1),
+                               (8,1,1),(8,7,1),(1,2,1),(2,3,1),
+                               (7,6,1),(6,5,1),(5,4,1),(3,4,1)])
+#2d graph Bridge-Top 
+c_tb = nx.Graph()
+c_tb.add_node(0,typ = "c1")
+c_tb.add_node(1,typ ="c2")
+c_tb.add_node(18,typ ="c2")
+c_tb.add_node(19,typ ="c2")
+c_tb.add_node(20,typ ="c2")
+c_tb.add_node(2,typ ="e1")
+c_tb.add_node(4,typ ="e1")
+c_tb.add_node(6,typ ="e1")
+c_tb.add_node(8,typ ="e1")
+c_tb.add_node(9,typ ="e2")
+c_tb.add_node(10,typ ="e1")
+c_tb.add_node(11,typ ="e2")
+c_tb.add_node(12,typ ="e1")
+c_tb.add_node(13,typ ="e2")
+c_tb.add_node(14,typ ="e1")
+c_tb.add_node(15,typ ="e2")
+c_tb.add_node(16,typ ="e1")
+c_tb.add_node(17,typ ="e2")
+c_tb.add_node(3,typ = "e2")
+c_tb.add_node(7,typ = "e2")
+for i in range(0,c_tb.number_of_nodes()):
+    c_tb.add_weighted_edges_from([(i,i+1,0)]) #no adjacent top-bridge
+    if i%2 == 0 and i<15:
+        c_tb.add_weighted_edges_from([(i,i+2,1)])#top-top bond type
+c_tb.add_weighted_edges_from([(0,i,1) for i in range(2,16,4)])#remaining top-top
+c_tb.add_weighted_edges_from([(0,l,0) for l in range(18,21)]) #Remaining adjacent sites
+c_tb.add_weighted_edges_from([(1,3,0),(3,5,0),(5,20,0),(20,1,0),(20,7,0),(7,9,0),(9,19,0),
+                             (11,19,0),(19,18,0),(18,13,0),(11,13,0),(1,18,0),(18,15,0),(15,17,0),(17,1,0)]) #no adjacent bridges
+c_tb.add_weighted_edges_from([(6,20,0),(18,14,0),(10,19,0)]) #remaining adjacent top-bridge
+
diff --git a/nmf/nmf.py b/nmf/nmf.py
index 0d907a3..c1127f7 100644
--- a/nmf/nmf.py
+++ b/nmf/nmf.py
@@ -55,7 +55,7 @@ def __init__(self,graph,num_ads,ads_names=None,ads_stoich=None,central_site=None
         self.make_a() #Energy of end atoms
         self.pprint() #Print details of class
         self.make_partition_function() 
-        #self.make_consistency_equations()
+        self.make_consistency_equations()
         
     def __repr__(self):
         return 'Cluster'
@@ -247,10 +247,16 @@ def make_consistency_equations(self):
         index = [self.names[i] for i in range(self.num_ads)]
         self.theta = pd.DataFrame(index=index,columns=self.end_sites)
         self.ln = 0
-        return
+        return self.const_eqns
     
     def get_theta(self, params, initial_guess=None, maxsteps=50, check_params=True, verify=False):
-
+        
+        if check_params:
+            key_list = list(params.keys())
+            for i in self.activities:
+                print(i)
+                assert str(i) in key_list, "Missing Input Values"
+      
         sub_eqns_expr = [None] * len(self.const_eqns)
         initial_guess1 = [0.5] * len(self.unknown_list)
         if initial_guess: