API tools faq
paste
Advertisement
shihab2196

C-H_Bond_Breaking_Location_PAO

shihab2196
Mar 26th, 2024
69
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.62 KB | None | 0 0
raw download clone embed print report
  1. # -*- coding: utf-8 -*-
  2. """
  3. Author: Shihab Ahmed
  4. Created on Mon Mar 25 14:26:26 2024
  5. """
  6.  
  7. import magnolia.bondfile_parser as bfp
  8. import numpy as np
  9. import pandas as pd
  10. import matplotlib.pyplot as plt
  11. import time
  12. import networkx as nx
  13. from networkx.algorithms import isomorphism
  14.  
  15. dirr      = r"C:\Users\arup2\OneDrive - University of California Merced\Desktop\LAMMPS\Base_Oil\PAO4\25_PAO4_200_O2_Soria\Production\1600"
  16. sim_dirr  = ['Sim-1', 'Sim-2', 'Sim-3']
  17.  
  18. bonddata = {}
  19. bo_cutoff = 0.3
  20. for sim in sim_dirr:
  21.     bondfilepath = dirr+'\\'+sim+'\\bonds.reaxc'
  22.     bonddata[sim]=bfp.parsebondfile(bondfilepath, cutoff=bo_cutoff,
  23.                                     bo=True, mtypes=True)
  24. #%%
  25.  
  26. ## Collect specific molecule_graph
  27. ## Here I am only collecting 75 PAO from three sim, 25 from each
  28. start = time.time()
  29.  
  30. mol_graphs  = []
  31. ref = None
  32.  
  33. for sim in sim_dirr:
  34.     bondfilepath = dirr+'\\'+sim+'\\bonds.reaxc'
  35.     fsbd     = bfp.parsebondfile(bondfilepath,firststep=True) #only the first step
  36.     fsneigh  = fsbd['neighbours'] # fsbd=first step bond data
  37.     atypes   = fsbd['atypes']    
  38.     fs_molecules= bfp.get_molecules(fsneigh)
  39.    
  40.     # get list of molecule_graph from the first steps
  41.     molecule_length = 92
  42.     fsneigh_graph = nx.Graph(fsneigh)
  43.     for molecule in fs_molecules:
  44.         if len(molecule)==molecule_length:
  45.             mol_graph = fsneigh_graph.subgraph(molecule)
  46.             mol_graphs.append(mol_graph)
  47.     break
  48.            
  49.  
  50. atom_map = {}
  51. site_map = {}
  52. if ref is None:
  53.     ref=mol_graphs[0]
  54.  
  55. for i, mol_graph in enumerate(mol_graphs):
  56.     print(i)
  57.     matcher = isomorphism.GraphMatcher(mol_graph, ref)
  58.     if matcher.is_isomorphic():
  59.         mapping = matcher.mapping
  60.        
  61.         for node in mol_graph.nodes:
  62.             atom_map[node] = mapping[node]
  63.            
  64.             # site mapping
  65.             if atypes[node] == 1: # H
  66.                 connected_carbon = list(mol_graph.neighbors(node))[0]
  67.                 site_map[node] = mapping[connected_carbon]
  68.             else:
  69.                 site_map[node] = mapping[node]
  70.  
  71. end = time.time()
  72. print(end-start)
  73. #%%------
  74. site_time = {}
  75. counter = 0
  76. data = []
  77. m_sim_type_check = []
  78. for sim in sim_dirr:
  79.     neighbors = bonddata[sim]['neighbours']
  80.     atypes    = bonddata[sim]['atypes']
  81.     mtypes    = bonddata[sim]['mtypes']
  82.    
  83.     seek = 'H61C30'
  84.     for step, neigh in neighbors.items():
  85.         molecules = bfp.get_molecules(neigh)
  86.         for molecule in molecules:
  87.             species = bfp.get_molecular_formula(molecule, atypes, 'HCO')
  88.             m_sim_type = sim+'-'+str(mtypes[list(molecule)[0]])
  89.             if species==seek and m_sim_type not in m_sim_type_check:
  90.                 data.append((molecule,step))
  91.                 m_sim_type_check.append(m_sim_type)
  92. #%%
  93. first = []
  94. for molecule,step in data:
  95.     fsbd     = bfp.parsebondfile(bondfilepath,firststep=True) #only the first step
  96.     fsneigh  = fsbd['neighbours'] # fsbd=first step bond data
  97.     atypes   = fsbd['atypes']    
  98.     fs_molecules= bfp.get_molecules(fsneigh)
  99.    
  100.     # get list of molecule_graph from the first steps
  101.     molecule_length = 92
  102.     fsneigh_graph = nx.Graph(fsneigh)
  103.     PAOs = []
  104.     for fs_molecule in fs_molecules:
  105.         if len(fs_molecule)==molecule_length:
  106.             PAOs.append(fs_molecule)
  107.    
  108.     for PAO in PAOs:
  109.         uncommon = PAO^molecule
  110.         if len(uncommon)==1:
  111.             single = list(uncommon)[0]
  112.             first.append(site_map[single])
  113.  
  114. plt.rc('font', size=14) # Default text sizes
  115. plt.rc('axes', titlesize=10) # Axes title size
  116. plt.rc('axes', labelsize=18) # Axes label size
  117. plt.rc('xtick', labelsize=16) # X-axis tick label size
  118. plt.rc('ytick', labelsize=16) # Y-axis tick label size
  119. plt.rc('legend', fontsize=12) # Legend fontsize
  120.  
  121. first = np.array(first)
  122. unique_elements, counts = np.unique(first, return_counts=True)
  123.  
  124. fig, ax = plt.subplots()
  125. bars = ax.bar(unique_elements-min(unique_elements)+1,counts,
  126.         color='tab:blue',edgecolor='black')
  127. ax.set_xlabel('H site')
  128. ax.set_ylabel('Count')
  129. ax.set_xticks(range(0,31,5))
  130. ax.set_ylim(top=6)
  131. ax.set_yticks(range(6))
  132.  
  133. # Adding a star sign on top of the highest bars
  134. max_height = max(counts)
  135. for bar in bars:
  136.     height = bar.get_height()
  137.     if height == max_height:
  138.         ax.text(bar.get_x() + bar.get_width() / 2., height,
  139.                 '$\star$', ha='center', va='bottom', fontsize=30, color='red')
  140.        
  141. plt.savefig(r'C:\Users\arup2\OneDrive - University of California Merced\Desktop\LAMMPS\Post Processing\lammps_processing\python_outputs\figures\H_breaking', dpi=300,
  142.             bbox_inches='tight')
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!