diff --git a/nmf/graphs.py b/nmf/graphs.py index 6a3034b..c6ca0e3 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 ="e1") 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,17 +44,64 @@ #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), (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 5524c8f..c1127f7 100644 --- a/nmf/nmf.py +++ b/nmf/nmf.py @@ -3,26 +3,51 @@ 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 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.end_sites = set([k['typ'] for i,k in self.graph.nodes(data=True)]) - if central_site: - self.end_sites.remove(central_site) + + #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.end_sites.remove("c") + 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)]) + 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 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 = {} + 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[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 @@ -36,12 +61,17 @@ 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=[] - 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 + ''' + 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 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 return def make_h(self): @@ -65,23 +95,29 @@ 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): ''' - Return list of dictionaries with length of list equal to number of adsorbates. - Each dictionary consists of symbol for each end node type. + 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=[] - 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) + self.a = {} + self.unknown_list = [] + 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'{etype}'] = df return def get_edge_value(self,u,v,w,G): @@ -103,25 +139,75 @@ 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 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 u,v,w in G.edges(data=True): - th *= self.get_edge_value(u,v,w,G) + for i,ads in enumerate(self.names): #loop over number of adsorbates + 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: + 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: + th=0 + else: + th=0 + + 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: + 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}") + print(f"z={self.total_expr}") return self.total_expr def make_consistency_equations(self): @@ -150,7 +236,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]]) @@ -161,15 +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): + 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: + for i in self.activities: + print(i) assert str(i) in key_list, "Missing Input Values" - - #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: @@ -208,24 +295,32 @@ 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("\nAds\tAds\tBond\tInteraction-Symbol") + 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}\t{self.ads_stoich[self.names.index(ads)]}") + + 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}") - 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\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): return cluster_draw(self.graph) + 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..0de37b5 100644 --- a/nmf/util.py +++ b/nmf/util.py @@ -89,11 +89,22 @@ 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) 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")