How to use the output from KSHELL

Load and view data from KSHELL

After running KSHELL, your work directory will look similar to this:

Ne20_usda.sh
Ne20_usda_m0p.wav
Ne20_usda_p.ptn
count_dim.py
kshell.exe
log_Ne20_usda_m0p.txt
log_Ne20_usda_tr_m0p_m0p.txt
save_input_ui.txt
transit.exe
usda.snt

All the level data are located in log_Ne20_usda_m0p.txt and all the transition data are located in log_Ne20_usda_tr_m0p_m0p.txt.

The log files are easily read with the kshell-utilities package. See the docstrings in the kshell-utilities repository for extended documentation. Install the package with pip:

pip install kshell-utilities

Create a blank Python file with your favourite editor. Lets name it ne20.py and lets place it in the results directory of the 20Ne calculation which is ~/kshell_results/ne20 according to this guide. However, the use of ~/ as a shortcut to your home directory is not standard in Python and is discouraged to be used, so if you wish to specify the path to your home directory, use the actual path. For macOS: /Users/<your username>/kshell_results/ne20. For most Linux distros: /home/<your username>/kshell_results/ne20. We use the loadtxt function to read the results from KSHELL:

import kshell_utilities as ksutil

def main():
    ne20 = ksutil.loadtxt(path=".")

if __name__ == "__main__":
    main()

The use of a name guard (if __name__ == "__main__":) is required because kshell-utilities uses Python's multiprocessing module which requires this to function properly. Note that path is a period (.). This simply means that the KSHELL results are located in the same directory as the Python file ne20.py. If we do not place ne20.py in the same directory as the KSHELL results, then we need to specify either the relative path to the KSHELL results from ne20.py or the absolute path of the KSHELL results which is (macOS) /Users/<your username>/kshell_results/ne20. Back to the loadtxt function. ne20 is an instance containing several useful attributes. To see the available attributes:

> print(ne20.help)
['debug',
'fname_ptn',
'fname_summary',
'gamma_strength_function_average_plot',
'gsf',
'help',
'level_density_plot',
'level_plot',
'levels',
'model_space',
'negative_spin_counts',
'neutron_partition',
'nucleus',
'parameters',
'path',
'proton_partition',
'transitions_BE1',
'transitions_BE2',
'transitions_BM1',
'truncation']

To see the energy, 2*angular momentum and parity of each level:

> print(ne20.levels)
[[-40.467   0.      1.   ]
[-38.771   4.      1.   ]
[-36.376   8.      1.   ]
[-33.919   0.      1.   ]
[-32.882   4.      1.   ]
[-32.107  12.      1.   ]
...
[-25.978  12.      1.   ]
[-25.904  10.      1.   ]
[-25.834   8.      1.   ]
[-25.829   2.      1.   ]]

Slice the array to get only selected values, if needed (ne20.levels[:, 0] for only the energies). To see 2*spin_initial, parity_initial, Ex_initial, 2*spin_final, parity_final, Ex_final, E_gamma, B(.., i->f), B(.., f<-i)] for the M1 transitions:

> print(ne20.transitions_BM1)
[[4.0000e+00 1.0000e+00 1.6960e+00 ... 7.5850e+00 5.8890e+00 0.0000e+00]
[4.0000e+00 1.0000e+00 1.6960e+00 ... 9.9770e+00 8.2810e+00 4.8200e-01]
[4.0000e+00 1.0000e+00 7.5850e+00 ... 9.9770e+00 2.3920e+00 1.1040e+00]
...
[4.0000e+00 1.0000e+00 1.3971e+01 ... 1.4638e+01 6.6700e-01 6.0000e-03]
[0.0000e+00 1.0000e+00 1.4126e+01 ... 1.4638e+01 5.1200e-01 2.0000e-02]
[2.0000e+00 1.0000e+00 1.4336e+01 ... 1.4638e+01 3.0200e-01 0.0000e+00]]

Visualise data from KSHELL

Create a level density plot

You can easily create a level density plot by

ne20.level_density_plot()

An alternative way is:

ground_state_energy: float = ne20.levels[0, 0]
ksutil.level_density(
    levels = ne20.levels[:, 0] - ground_state_energy,
    bin_width = 0.2,
    plot = True
)

Note that scaling the excitation energies by the ground state energy is required with this method. If you want greater control of matplotlib plotting parameters, use this method:

import matplotlib.pyplot as plt

ground_state_energy: float = ne20.levels[0, 0]
bins, density = ksutil.level_density(
    levels = ne20.levels[:, 0] - ground_state_energy,
    bin_width = 0.2,
    plot = False,
)
plt.step(bins, density)
plt.show()

The bin_width is in the same energy units as your results, which for KSHELL is MeV. The two latter ways of generating the plot does not require that the data comes from KSHELL. Use any energy level data normalised to the ground state energy. The plot will look like this:

Click to see level density plot

Create a level plot / level scheme

To generate a level plot:

ne20.level_plot()

or

import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ksutil.level_plot(
    levels = ne20.levels,
    ax = ax
)
plt.show()

Click to see level plot

Both ways of generating the level plot supports selecting what total angular momenta to include in the plot, and how many levels per angular momentum.

ne20.level_plot(
    include_n_levels = 3,
    filter_spins = [0, 3, 5]
)

Click to see filtered level plot

Create a gamma strength function plot

The gamma strengh function (averaged over total angular momenta and parities) can easily be calculated in several ways. The quickest way is

ne20.gsf()

which is an alias for the following function call:

ne20.gamma_strength_function_average_plot(
    bin_width = 0.2,
    Ex_max = 5,
    Ex_min = 20,
    multipole_type = "M1",
    plot = True,
    save_plot = False
)

The default parameters are applied if no function arguments are supplied. If you want to have greater control over the plotting procedure, then this solution is better:

import matplotlib.pyplot as plt

bins, gsf = ne20.gamma_strength_function_average_plot(
    bin_width = 0.2,
    Ex_max = 50,
    Ex_min = 5,
    multipole_type = "M1",
    plot = False,
    save_plot = False
)
plt.plot(bins, gsf)
plt.show()

since you yourself have control over the matplotlib calls. Note that Ex_max is set to higher energy than what you usually get from the KSHELL calculations. The default upper limit is set large as to include all levels of any KSHELL calculation. The final way of doing it is:

import matplotlib.pyplot as plt

bins, gsf = ksutil.gamma_strength_function_average(
    levels = ne20.levels,
    transitions = ne20.transitions_BM1,
    bin_width = 0.2,
    Ex_min = 5,
    Ex_max = 20,
    multipole_type = "M1"
)
plt.plot(bins, gsf)
plt.show()

where the difference is that you supply the levels and transitions arrays. I'd not recommend this final solution unless you have level and transition data from some other place than KSHELL. The parameters bin_width, Ex_max and Ex_min are in the same unit as the input energy levels, which from KSHELL is in MeV. bin_width is the width of the bins when the level density is calculated. Ex_min and Ex_max are the lower and upper limits for the excitation energy of the initial state of the transitions.

Click to see gamma strength function plot