Skip to content

Commit

Permalink
scripts sync
Browse files Browse the repository at this point in the history
  • Loading branch information
Phoelionix committed Jan 18, 2025
1 parent ca0a861 commit a4d5863
Show file tree
Hide file tree
Showing 13 changed files with 1,084 additions and 205 deletions.
6 changes: 3 additions & 3 deletions config/version.txt
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
Version 2: T.N.T.
Totally Non-Thermal electron plasma simulation
Alexander Koszlov and Alaric Sanders 2020
Version 3: [Acronym]
[Name]
Alexander Koszlov, Alaric Sanders and Spencer Passmore 2024
37 changes: 15 additions & 22 deletions input/debug/Debug.mol
Original file line number Diff line number Diff line change
@@ -1,37 +1,38 @@
#ATOMS
C 990
C 615.16
N 195.16
O 887.2
S 10
Gd_fast 3.08


#BOUND_FREE_EXCLUSIONS
S

#VOLUME
40000 // Volume per molecule in Angstrom^3.
30820 // Volume per molecule in Angstrom^3.
500 // Length scale of sample in Angstrom. Used for effective escape rate of photo-electrons (assumes surrounded by a void).
none // Spatial boundary shape - options are none, spherical, cylindrical, planar


#PULSE
7112 // Photon energy in eV.
15 // Pulse width in femtoseconds (defined as FWHM for Gaussian pulse).
gaussian // Pulse shape

#PROBE_PULSE
25

#USE_COUNT
true // active (using count)?
0.1 // Photon density (x10^12 ph.µm-2) (3.5e12 um-2)
1.75 // Photon density (x10^12 ph.µm-2) (3.5e12 um-2)

#NUMERICAL
2000 // Initial guess for number of time step points for rate equation.
21 // Number of threads in OpenMP.
21 // Number of threads in OpenMP.

#DYNAMIC_GRID
M // Grid regions preset, options are: dismal, low, medium, high, (referring to accuracy), and no_dirac (static region is used rather than dynamic dirac regions)
5.05 // Grid update period in fs, (dynamic grid only).
2 // Time before transforming initial guess grid to dynamic grid in fs.

#OUTPUT
2400 // Number of time steps in the output files.
800 // Number of time steps in the output files.
4000 // number of free-electron grid points in output file.
N // Write atomic charges in a separate file (Y/N)?
Y // Write intensity in a separate file (Y/N)?
Expand All @@ -40,20 +41,12 @@ Y // Write data for molecular dynamics (MD) in a separate file (Y/N)?
#DEBUG
99 //10 // Simulation cutoff time in fs
0.01 // Interval to update current timestep [fs].
0 // Interval **in steps** between live plot (_live_plot.png) updates
90 // Period between backups in minutes
####END#####

Notes:

Parameters approximate (pump) pulse used in Nass et al 2020 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7156470/
Hen egg-white lysozyme https://files.rcsb.org/pub/pdb/validation_reports/et/4et8/4et8_full_validation.pdf
Ignores:
- precipitant (NaCl)
- surrounding water, which would filter photoelectron distribution in crystal with own.
- If photoelectrons aren't slowed by boundary significantly, this is on a timescale shorter
than the pulse.
- spatial boundary, since the water droplet size is unknown

from many runs, distributions are very well-approximated by low accuracy grid,
for the purposes of both a) atomic form factor determination and b) debye length
determination for hybrid models.
see lys_Gd_salt_solvated

ignoring volume change from removing salt (removing salt probably increases the volume negligibly)
14 changes: 12 additions & 2 deletions run_sims.sh
Original file line number Diff line number Diff line change
Expand Up @@ -20,11 +20,21 @@ set -x

#############################

#./ac4dc_split input/lysozyme/lys_Gd_salt_solvated_fast
#./ac4dc_split input/lysozyme/lys_solvated_light
# ./ac4dc_split input/lysozyme/lys_Gd_solvated_fast

# ./ac4dc input/lysozyme/lys_solvated
# ./ac4dc input/lysozyme/lys_Gd_salt_solvated
# ./ac4dc input/lysozyme/lys_Gd_solvated

./ac4dc input/lysozyme/lys_Gd_solvated
./ac4dc input/lysozyme/lys_Gd_salt_solvated # if above works there might be a new memory leak in the code
#./ac4dc input/nass/nass_probe_0
#./ac4dc input/nass/nass_probe_35
#./ac4dc input/nass/nass_probe_62
./ac4dc input/I3C/I3C_5fs_olaf
./ac4dc input/I3C/I3C_20fs_olaf
# ./ac4dc input/I3C/I3C_5fs_olaf
# ./ac4dc input/I3C/I3C_20fs_olaf



Expand Down
84 changes: 57 additions & 27 deletions scripts/_generate_plots.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,8 +25,10 @@
####
ELECTRON_DENSITY = False # Whether to use electron density for free distribution plots. Energy density if False
###
PLOT_ELEMENT_CHARGE= False #
PLOT_FREE_CONTINUUM = True
PLOT_ELEMENT_CHARGE= True #
PLOT_FREE_CONTINUUM = False
PLOT_SPLIT_FREE_CONTINUUMS = False
PLOT_COMBINED_SPLIT_FREE_CONTINUUMS = False
PLOT_FREE_SLICES=False
PLOT_ION_RATIOS=False
PLOT_ION_RATIOS_BARS= False
Expand All @@ -39,12 +41,15 @@
# FIGWIDTH = COLUMNWIDTH#/2
# FIGHEIGHT = FIGWIDTH*1/2#*9/16

FIGWIDTH = COLUMNWIDTH*1.1
FIGHEIGHT = FIGWIDTH*9/16
#FIGWIDTH = COLUMNWIDTH*1.1
#FIGHEIGHT = FIGWIDTH*9/16
COLWIDTH = COLUMNWIDTH*1.1
ROWHEIGHT = COLWIDTH*9/16

# FIGWIDTH = COLUMNWIDTH/2
# FIGHEIGHT = FIGWIDTH*9/16

#DPI = 100
DPI = 800
##
END_T = None
Expand All @@ -70,15 +75,20 @@ def main():
assert valid_folder_names, "One or more arguments (directory names) were not present in the output folder."
for data_folder in sys.argv[1:]:
label = data_folder +'_Plt'
make_some_plots(data_folder,molecular_path,label,dname_Figures,PLOT_ELEMENT_CHARGE,PLOT_ION_RATIOS,PLOT_FREE_CONTINUUM,PLOT_FREE_SLICES,PLOT_ION_RATIOS_BARS,PLOT_ORBITAL_DENSITIES,PLOT_PHOTO_RATES)
make_some_plots(data_folder,molecular_path,label,dname_Figures,PLOT_ELEMENT_CHARGE,PLOT_ION_RATIOS,PLOT_FREE_CONTINUUM,PLOT_FREE_SLICES,PLOT_ION_RATIOS_BARS,PLOT_ORBITAL_DENSITIES,PLOT_PHOTO_RATES,PLOT_SPLIT_FREE_CONTINUUMS,PLOT_COMBINED_SPLIT_FREE_CONTINUUMS)

def make_some_plots(mol_name,sim_output_parent_dir, label,figure_output_dir, tot_charge=False,bound_ionisation=False,free=False,free_slices=False,bound_ionisation_bar=False,orbital_densities_bar=False,photo_rates = False):
def make_some_plots(mol_name,sim_output_parent_dir, label,figure_output_dir, tot_charge=False,bound_ionisation=False,free=False,free_slices=False,bound_ionisation_bar=False,orbital_densities_bar=False,photo_rates = False,split_free=False,combined_split_free=False):
'''
Arguments:
mol_name: The name of the folder containing the simulation's data (the csv files). (By default this is the stem of the mol file.)
sim_data_parent_dir: absolute path to the folder containing the folders specified by target_handles.
'''

# extra homeless options
load_specific_atoms = ["C","N","O","S","Gd_fast"]#["Gd_fast"] #["C"]#["Fe_singleShell","C"] #None #["C","N","O"] #If plotting free dsitribution, will plot only those specified here.
split_continuums_to_load = [] #["Gd_fast"] #["C"]#["Fe_singleShell","C"] # None is not valid. Use "[]"


############
# File/directory names
#######
Expand All @@ -87,25 +97,42 @@ def make_some_plots(mol_name,sim_output_parent_dir, label,figure_output_dir, tot
fname_free = "free"
fname_HR_style = "HR_style"
fname_bound_dynamics = "bound_dynamics"
load_specific_atoms = ["C"] #["C"]#["Fe_singleShell","C"] #None #["C","N","O"] #If plotting free dsitribution, will plot only those specified here. If None is specified, will just use the full continuum freeDist.csv.




##############
if combined_split_free or split_free:
assert split_continuums_to_load == "all" or len(split_continuums_to_load) > 0

if load_specific_atoms is not None:
label+="_"
for elem in load_specific_atoms:
label+=elem
pl = Plotter(mol_name,sim_output_parent_dir,use_electron_density = ELECTRON_DENSITY,end_t = END_T,
initialise=False)
pl.get_atoms(load_specific_atoms,split_continuums_to_load) # so that plotter knows num continuums


pl = Plotter(mol_name,sim_output_parent_dir,use_electron_density = ELECTRON_DENSITY,end_t = END_T,load_specific_atoms=load_specific_atoms)
num_atoms = len(pl.statedict)
num_subplots = tot_charge + free_slices + bound_ionisation_bar + (bound_ionisation+ orbital_densities_bar+ photo_rates)*num_atoms
num_subplots+= free*pl.num_continuums()
#num_subplots+= free
pl.setup_axes(num_subplots)
num_subplots = tot_charge + free_slices + bound_ionisation_bar + (bound_ionisation+ orbital_densities_bar+ photo_rates)*num_atoms + free + combined_split_free + split_free*pl.num_continuums()


if free:
pl.initialise(None,num_subplots,"full") # Load full continuum.
pl.plot_free(log=True,ylog=False,cmin=10**(-8),cmax=10**(-3.609),ylim=[0,8000],keV=True,show_title=False)
pl.update_inputs(load_specific_atoms=load_specific_atoms,split_continuums_to_load=split_continuums_to_load)
else:
pl.initialise(load_specific_atoms,num_subplots,split_continuums_to_load,)

if num_subplots > 1:
pl.fig.tight_layout()
pl.fig.subplots_adjust(left=0.12/pl.axs.shape[0], bottom=None, right=None, top=None, wspace=0.2, hspace=None)
pl.fig.subplots_adjust(left=0.12/pl.axs.shape[0], bottom=None, right=None, top=None, wspace=0.2, hspace=0.2)#hspace=None)

if tot_charge:
#NOTE ensure load_specific_atoms is None or does not exclude `atoms` if `atoms` is passed.
pl.plot_tot_charge(ylim=[None,None],every=1,charge_difference=False,legend_loc="best")
#pl.plot_tot_charge(ylim=[None,None],xlim=[-18,18],every=1,charge_difference=False,legend_loc="best")
#pl.plot_tot_charge(ylim=[None,None],every=1,charge_difference=True,legend_loc="best",atoms=["C","N","O"])
#pl.plot_tot_charge(ylim=[0,6],every=1,charge_difference=False,legend_loc="best",atoms=["C","N","O"])
#pl.plot_tot_charge(ylim=[0,6],plot_legend=False,every=1,charge_difference=True,legend_loc="best")
Expand Down Expand Up @@ -142,17 +169,25 @@ def make_some_plots(mol_name,sim_output_parent_dir, label,figure_output_dir, tot
# pl.plot_charges_leonov_style("Si",show_pulse_profile=True,xlim=[-20,20],ylim=[0,1.099])
# pl.fig.set_figwidth(6.662*0.7)
# pl.fig.set_figheight(6*0.7)

if free:
every_e = 1 # every nth energy plotted.
every_t = 1

if combined_split_free:
pl.plot_free(log=True,ylog=False,cmin=10**(-8),cmax=10**(-3.609),ylim=[0,8000],keV=True,show_title=False,every=every_t,every_e=every_e)
if split_free:
#pl.plot_free(log=True,cmin=10**(-7.609),cmax=1e-3,ylim=[10,8000])
for _c in range(pl.num_continuums()):
pl.plot_free(log=True,ylog=False,cmin=10**(-8),cmax=10**(-3.609),ylim=[0,8000],keV=True,continuum=_c)
pl.fig.tight_layout()
for _c in range(pl.num_continuums()-2):
print(_c)
pl.plot_free(log=True,ylog=False,cmin=10**(-8),cmax=10**(-3.609),ylim=[0,8000],keV=True,continuum=_c,every=every_t,every_e=every_e)



# #Leonov
# ymax = 9e3
# pl.plot_free(log=True, cmin=10**(-6.609),cmax = 10**(-2), every=5,mask_below_min=True,cmap='turbo',ymax=ymax,leonov_style=True)
# pl.fig.set_figwidth(6.662*0.7*1.16548042705)
# pl.fig.set_figheight(6*0.7)

if free_slices:
pl.initialise_step_slices_ax()
from plotter_core import fit_maxwell, maxwell
Expand Down Expand Up @@ -290,17 +325,12 @@ def make_some_plots(mol_name,sim_output_parent_dir, label,figure_output_dir, tot
# pl.ax_steps.set_xlim([0,1800])
# pl.ax_steps.set_ylim([0.5e-4,0.5])
pl.delete_remaining_axes()
plt.gcf().set_figheight(FIGHEIGHT*pl.axs.shape[1])
plt.gcf().set_figwidth(COLWIDTH*pl.axs.shape[0])
plt.gcf().set_figheight(ROWHEIGHT*pl.axs.shape[1])


if num_subplots > 2:
plt.gcf().set_figwidth(FIGWIDTH*pl.axs.shape[0])
# elif num_subplots == 2:
# plt.gcf().set_figheight(FIGHEIGHT*pl.axs.shape[0])
else:
plt.gcf().set_figwidth(FIGWIDTH)
plt.gcf().set_figheight(FIGHEIGHT)
#plt.tight_layout()
plt.savefig(figure_output_dir + label + figures_ext,dpi=DPI)
plt.savefig(figure_output_dir + label + figures_ext,dpi=DPI,bbox_inches='tight')
plt.close()

if __name__ == "__main__":
Expand Down
117 changes: 117 additions & 0 deletions scripts/element_ion_plot.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@
# Comparing all elements between two simulations.

import matplotlib
matplotlib.use("pgf")
matplotlib.rcParams.update({
"pgf.texsystem": "pdflatex",
'font.family': 'serif',
'font.size': 10,
'text.usetex': True,
'pgf.rcfonts': False,
})
import matplotlib.pyplot as plt
import numpy as np
from plotter_core import Plotter
import sys, traceback
import os.path as path
import os
from QoL import set_highlighted_excepthook

OCCUPANCY = True # ONLY IMPLEMENTED FOR PLOT MODE 1 CURRENTLY
CHARGE_DIFFERENCE = False # Set True if want initially ionised species' traces to start from origin
PLOT_DERIVATIVE = False # Plot the rate of avg charge gain
#YLIM = [0,8]
PLOT_MODE = 1 # 1: plot element total charges # 2: plot orbital charges
YLIM = [0,None]
SCALE = 6
FIGWIDTH = SCALE*2/3
FIGHEIGHT = SCALE/2
#XLIM = [None,None]
XLIM = [None,None]
#YLIM=[0,6]
#YLIM=[20,24]

ATOM = "C"
#ATOMS = ("Cr_LDA",)

CUSTOM_LEGEND = None
#CUSTOM_LEGEND=("1s: CNO","2s: CNO","2p: CNO","1s: CNO,S,Gd","2s: CNO,S,Gd","2p: CNO,S,Gd")
#CUSTOM_LEGEND = ("Primary ionization only", "All ionization",)
def main():
set_highlighted_excepthook()


# Basic num arguments check
assert len(sys.argv[1:]) > 0, "Usage: python compare_ion.py <sim_handle_1> <sim_handle_2> ..."

molecular_path = path.abspath(path.join(__file__ ,"../../output/__Molecular/")) + "/"
dname_Figures = "../../output/_Graphs/plots/"
dname_Figures = path.abspath(path.join(__file__ ,dname_Figures)) + "/"
valid_folder_names= True
for i, data_folder in enumerate(sys.argv[1:]):
if not path.isdir(molecular_path+data_folder):
valid_folder_names = False
print("\033[91mInput error\033[0m (argument \033[91m"+str(i)+ "\033[0m): folder name not found.")
assert valid_folder_names, "One or more arguments (directory names) were not present in the output folder."
data_folders = sys.argv[1:]
label = data_folders[0] + "_"+data_folders[1]
make_some_plots(data_folders,molecular_path,label,dname_Figures,plot_derivative=PLOT_DERIVATIVE)

def make_some_plots(mol_names,sim_output_parent_dir, label,figure_output_dir,plot_derivative = False):
'''
Arguments:
mol_name: The name of the folder containing the simulation's data (the csv files). (By default this is the stem of the mol file.)
sim_data_parent_dir: absolute path to the folder containing the folders specified by target_handles.
'''

############
# File/directory names
#######
figures_ext = "" #.png
for plot_mode in (PLOT_MODE,):
fig, axs = plt.subplots(3, 3, sharey=True, facecolor='w')



cmap = plt.get_cmap("Dark2")

# Hacky way to get overlaid plots TODO
pl = Plotter(mol_names[0],sim_output_parent_dir)
pl.setup_axes(1)
for m, mol_name in enumerate(mol_names):
pl.__init__(mol_name,sim_output_parent_dir)
pl.num_plotted = 0 # ƪ(˘へ˘ ƪ)
if plot_mode == 1:
colours = [cmap(m)]
pl.plot_tot_charge(every=1,colours = colours,atoms = [ATOM],plot_legend=(m==0),xlim=XLIM,ylim=YLIM,charge_difference=CHARGE_DIFFERENCE,plot_derivative=PLOT_DERIVATIVE,occupancy=OCCUPANCY)

if plot_mode == 2:
ax = pl.plot_orbitals_charge(every=1,atom = ATOM,plot_legend=False,xlim=XLIM,ylim=YLIM,plot_derivative=PLOT_DERIVATIVE)


ax = pl.axs[0][0]
if CUSTOM_LEGEND is None:
ax.legend(bbox_to_anchor=(1.02, 1),loc='upper left', ncol=1,handlelength=1) # Top right legend.
else:
handles,_ = ax.get_legend_handles_labels()
handles = list(handles)
assert len(CUSTOM_LEGEND)==len(handles)
ax.legend(handles,CUSTOM_LEGEND,bbox_to_anchor=(1.02, 1),loc='upper left', ncol=1,handlelength=1)
plt.gcf().set_figwidth(FIGWIDTH)
plt.gcf().set_figheight(FIGHEIGHT)
#plt.gcf().tight_layout()
#plt.tight_layout()

if plot_mode == 1:
qualifier = f"_{ATOM}-Comp"
if plot_mode == 2:
qualifier = f"_{ATOM}-OrbsComp"
if PLOT_DERIVATIVE:
qualifier+="-deriv"
plt.savefig(figure_output_dir + label +qualifier + figures_ext,bbox_inches='tight')
plt.close()

if __name__ == "__main__":
main()

#TODO: Change size to match my screen by default, add --y option
2 changes: 1 addition & 1 deletion scripts/generate_interactive.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,7 @@
P.ELECTRON_DENSITY = False # if False, use energy density
P.ALSO_MAKE_PLOTS = False # Generate static plots for each simulation using _generate_plots.py
P.SINGLE_FRAME = False # Save a png using plotly .
P.NAMING_MODE = 0 # For legend. 0: full details of sim parameters + sim name | 1: elements in sim|
P.NAMING_MODE = 1 # For legend. 0: full details of sim parameters + sim name | 1: elements in sim|
P.SCALE_DENSITY_BY_THOUSAND = False # Use cubic nm rather than cubic angstrom for measuring energy density
P.END_T = 9999 # Put at value to cutoff times early. Note if multiple handles inputted will cutoff all to earliest end time.
P.POINTS = 70
Expand Down
Loading

0 comments on commit a4d5863

Please sign in to comment.