API tools faq
paste
Advertisement
shihab2196

heatmap_species_molecule_B

shihab2196
Jan 11th, 2024
93
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 4.27 KB | None | 0 0
raw download clone embed print report
  1. # -*- coding: utf-8 -*-
  2. """
  3. Author: Shihab Ahmed
  4. Created on Mon Oct 23 05:15:03 2023
  5. """
  6.  
  7. ############################################
  8. ##
  9. #   Heatmap of most abundance species
  10. #   Abundance based on mean number of species
  11. #   at the last 100(or anything you want) ps.
  12. #   Basically, it show only the stable products not intermediates
  13. #   Using species file
  14. ##
  15. ############################################
  16.  
  17. import magnolia.bondfile_parser as bfp
  18. import matplotlib.pyplot as plt
  19. from matplotlib.colors import LogNorm
  20. import random
  21. import pandas as pd
  22. import magnolia.speciesfile_parser as sfp
  23. import seaborn as sns
  24. import numpy as np
  25.  
  26.  
  27. # User Inputs
  28. savedir = r'C:\Users\arup2\OneDrive - University of California Merced\Desktop\LAMMPS\Post Processing\lammps_processing\python_outputs\heatmap_of_species'
  29.  
  30. commmon = r"C:\Users\arup2\OneDrive - University of California Merced\Desktop\LAMMPS\Antioxidants\ABCDE\B\B_300_O2\Production\1250"
  31.  
  32. ### Very Important !!!!!!!!!!!!!!!!!!!!!!!!!!
  33. # if isSpecefile==True: use species file, else use bond file
  34. isSpeciesFile = True
  35.  
  36. bo_cutoff = 0.30
  37. flag = True
  38. simulations = ['\\Sim-1']#,'\\Sim-2','\\Sim-3']
  39. for sim in simulations:
  40.     directory = commmon+sim
  41.     speciesfile_path = directory+'\\species.out' if bo_cutoff==0.3 else directory+f'\\species_{bo_cutoff}.out'
  42.     bondfile_path    = directory+'\\bonds.reaxc'
  43.    
  44.     if isSpeciesFile:
  45.         # using species file (too fast)
  46.         current = sfp.get_species_count(speciesfile_path)
  47.     else:
  48.         # using bond file (too slow)
  49.         current = bfp.get_species_count(bondfile_path,
  50.                                         atomsymbols=['H','C','O'],
  51.                                         cutoff=bo_cutoff)
  52.    
  53.     if flag:
  54.         output_species = current
  55.         flag = False
  56.     else:
  57.         output_species = output_species.add(current,fill_value=0)
  58.  
  59. ## other wise plot doesn't end at 1000 ps
  60. if output_species.columns[0]!= 0:
  61.     output_species.insert(0, 0, output_species[output_species.columns[0]])
  62. #%%
  63. # User Inputs
  64. timestep = 0.25
  65. skipts   = 0#(1150-300)/4
  66. nspecies = 10
  67. tailps   = 10  #ps
  68. lignin   = {"A":'H30C19O3',"B":'H34C26O4'}
  69. exclude  = ['O2',lignin['B']]
  70. n_xticks = 4  # how many xtick division in x axis i.e. 1+n_xticks ticklabels
  71. nsim     = len(simulations)
  72. fontsize = 10
  73.  
  74. inputfile_name= "species.out" if isSpeciesFile else "bonds.reaxc"
  75. details  = ('\n\nThis heatmap is generated using {} file. Number of '
  76.             'species are average\nover {} simulations.'
  77.             'Species are sorted by mean species count of last {} ps.\n'
  78.             'bo_cutoff {} is used here, '
  79.             'First {} species are showing here.\n'.format(inputfile_name, nsim,
  80.                                                           tailps, bo_cutoff,
  81.                                                           nspecies))
  82. title    = commmon[commmon.find('ABCDE'):]+ details
  83.  
  84.  
  85.  
  86. species = output_species.copy()/nsim # copy is important
  87. # species.replace(0,1e-20,inplace=True) #to avoid log(0) if use log scale
  88. species.columns = species.columns*timestep/1000-skipts
  89. # species    = species.loc[:,0:] # skip the ramping steps
  90. species    = species.drop(exclude)
  91.  
  92. ## mean of last {tailps} ps
  93. means          = species.loc[:,species.columns.max()-tailps:].mean(axis=1)
  94. sorted_index   = means.sort_values(ascending=False).index.tolist()
  95. species        = species.loc[sorted_index]
  96.  
  97. species       = species.head(nspecies)
  98. species.index = bfp.make_molecular_formula_latex(species.index,sort=True)
  99.  
  100. xticks_freq = int(species.columns.size/n_xticks)
  101. norm = None#LogNorm(vmin=1, vmax=species.max().max()) #log-scale or set None
  102. fig, ax = plt.subplots(figsize=(10,5))
  103. sns.heatmap(species,cmap='jet',xticklabels=xticks_freq, ax=ax,
  104.             cbar_kws={'label': 'Number of species'},norm=norm)
  105.  
  106. # make the xticks integers
  107. xticklabels = ax.get_xticklabels()
  108. xticks_ = np.array([float(label.get_text()) for label in xticklabels])
  109. xticks  = xticks_.astype(int)
  110. ax.set_xticklabels(xticks,fontsize=fontsize)
  111. # ax.set_xlim(0,100*4)
  112.  
  113. ax.set_xlabel('Time (ps)',fontsize=fontsize)
  114. plt.tick_params(left=False)
  115. plt.title(title,fontsize=10)
  116. plt.yticks(rotation=0)
  117. plt.savefig(savedir+'\\hm_species_{}'.format(random.randint(0,9999999999)),dpi=1000, bbox_inches='tight')
  118. plt.show()
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!