Table of Contents
Welcome to The Tech Basement!
The Tech Basement is my personal knowledge base! It is organised as a wiki, but I am also using it like a blog. The subjects are programming, maths and physics, which are subjects highly related to my the PhD I'm taking in computational nuclear physics.
Check out the latest entries in The Tech Basement below!
[SOLVED] Possible bug in the OBTD calculations
So there is a possible bug that I am trying to wrap my head around. In the OBTD log files from KSHELL, the OBTDs are listed in blocks where each block represents an initial and a final state aka. one specific transition. For example
w.f. J1= 0/2( 3) J2= 2/2( 1) B(L;=>), B(L ;<=) 0.00622 0.00207 <||L||> 3 1 -0.13590 0.11892 i j OBTD <i||L||j> OBTD*<||> 1 1 0.00244 5.79655 0.01416 1 2 -0.00012 1.54919 -0.00018 1 9 0.00000 0.00000 0.00000 1 10 0.00000 0.00000 0.00000 ...
This block represents the OBTDs for the transition from the 3rd 0- state to the 1st 1- state (the parity is not listed in the block, but I know it is negative). Now, there is no guarantee that the 3rd 0- state has higher energy than the 1st 1- state. It might be so, but it is not guaranteed. If the initial state however has lower excitation energy than the final state then the OBTD indices in the blocks ($i, j$) have to be swapped in order to make it so that the final state is of lower energy.
So I made a short script which runs through all of the OBTD blocks and looks up the excitation energy of the initial and final state for the given block. The result was that 628337 of 1280000 (49.09 %) blocks had initial state excitation energy lower than the final state. Thats good, as expected, some of the final states are of lower energy and some are of higher. The puzzle however is that when I make a sub-set of transitions, selecting only the transitions which are present within the $E_\gamma = [0, 3]$ MeV region of the $M1$ GSF then only 3 of 230397 (0.00 %) transitions have initial energies lower than final. I don't understand where in the process of calculating the GSF that this kind of energy sorting happens.
Visualising OBTDs
The one-body transition density (OBTD) is used to calculate the transition probability from an initial to a final state:
$$ \langle \Psi_f | \hat{O}_{\lambda \mu} | \Psi_i \rangle = \sum_{\alpha \beta} \langle \alpha | \hat{o}_{\lambda \mu} | \beta \rangle \langle \Psi_f | \hat{c}^\dagger_\alpha \hat{c}_\beta | \Psi_i \rangle \qquad(0) $$
with the OBTD defined as:
$$ \rho_{fi}(\alpha, \beta) = \langle \Psi_f | \hat{c}^\dagger_\alpha \hat{c}_\beta | \Psi_i \rangle. \qquad(1) $$
Assume $\alpha$ and $\beta$ ($i$ and $j$ in the list below) to be indices of the orbitals (not the $m$ substates) in the model space, which for sd-pf-sdg are 12 proton orbitals and 12 neutron orbitals. Let's select a subset of transitions so that we have transitions of gamma energies only within a certain gamma energy interval $[E_\text{min}, E_\text{max}]$. This might for example be the energy interval of the LEE. Let's also select transitions based on their multipolarities so that we look at only, for example, $M1$ transitions. Transitions of $0^- \rightarrow 1^-$ are within these requirements, and the raw OBTD data for one such transition is listed below. Note that the $M1$ operator has an $L$ and an $S$ term
$$ \hat{M1} = g_l \hat{L} + g_s \hat{S}.\qquad(2) $$
The file below lists some $L$ matrix element values for vanadium 50 and there is a corresponding file with the $S$ values.
w.f. J1= 0/2( 3) J2= 2/2( 1) B(L;=>), B(L ;<=) 0.00622 0.00207 <||L||> 3 1 -0.13590 0.11892 i j OBTD <i||L||j> OBTD*<||> 1 1 0.00244 5.79655 0.01416 1 2 -0.00012 1.54919 -0.00018 1 9 0.00000 0.00000 0.00000 1 10 0.00000 0.00000 0.00000 1 11 0.00000 0.00000 0.00000 ...
Click here to see the complete L list
Click here to see the complete S list
Since there are no transitions between proton and neutron orbitals (gamma decay cannot turn protons into neutrons), we can organise the OBTDs into two 12×12 grids, one grid for protons and one grid for neutrons, for example like this (protons):
Notice that the OBTD can be negative. If we don't care about the sign itself, we might take the absolute value of the OBTDs to get a feel of how large the overlap of the initial and final states are given the single-particle transition from $\alpha$ ($i$) to $\beta$ ($j$). One such 12×12 proton OBTD grid can be made for not just this one transition, but one can be made for each of the transitions we selected, which might be tens or hundreds of thousands of transitions. Let's then consolidate all of the thousands of 12×12 grids into a single 12×12 grid by adding the grids at each grid value. This results in some very large grid values whose magnitude doesn't really matter; what matters is how large the values are relative to each other. Let's then normalise by dividing each grid value by the sum of all the grid values and multiply by 100 to express them as a percentage of the sum of all the OBTDs:
Note that the percentages do not add to 100% because they are normalised to both proton and neutron values and the above plot shows only proton values. The sum of the proton and neutron grid values add up to 100%.
An important consideration to make is that the OBTD eq. (1) has no dependency on any transition operator. The OBTD simply tells us the overlap between the initial and final wavefunctions given a single-particle transition in the initial wavefunction. The OBTD does not care about electromagnetic transition rules, so if we want to consider OBTDs for a subset of transitions with an $M1$ characteristic we need to do some filtering.
Note that the $M1$ operator in eq. (2) consists of an $L$ and an $S$ term, the values for $L$ and $S$ for a single transition are given in the above lists. We see that in some cases a non-zero OBTD value has corresponding $L$ and $S$ values which are zero. In such a case, I think it would be wrong to include that OBTD value in the 12×12 grid since its contribution is zeroed by transition rules. Here I'm assuming that the angular momentum and parity transition rules manifest as the $L$ and $S$ terms being zero. If you look at eq. (0), where else could the transition rules come into play? There are only two terms.
Now what can we say about these two latest 12×12 grids… First, the single-particle transitions seem to respect the transition rules for $M1$ transitions. Recall that $M1$ transitions have no change in parity from the initial to the final state which means that the only allowed single-particle transitions are within each major shell. We can clearly see block structures within $sd$, $pf$, and $sdg$. Second, $M1$ transitions allow only $\Delta j = 1$, recall the rule
$$ \Delta j = | j_i - j_f |, ..., j_i + j_f $$
which means that, for example, a single-particle transition from $s1/2$ to $d5/2$ has a minimum $\Delta j$ of 2, and consequently that this single-particle transition is illegal for $M1$ transitions.
We can see that even before explicitly respecting the $M1$ transition rules by skipping OBTDs where $L$ and $S$ are zero, the OBTDs in the 12×12 grid already agree with $\Delta j = 1$ and $\Delta \pi = \text{no}$. This probably comes from the fact that I already made a selection of transitions which are $M1$ transitions, but just because I selected $M1$ transitions does not mean that all single-particle transitions within the $M1$ transitions (eq. (0)) are allowed.
A still open question to me is why the $1s1/2 \rightarrow 0d3/2$ OBTDs disappeared when
$$ \Delta j = 3/2 - 1/2, ..., 3/2 + 1/2 = 1, 2 $$
allows $M1$ transitions…
Giant Resonances
This article contains my notes about the giant resonances in nuclear physics. The information is mainly from Giant Resonances: Fundamental High-Frequency Modes of Nuclear Excitation by M. N. Harakeh. Please see that book for more details!
Overview of effective interactions
In this article you'll find an overview of many effective interactions. My goal is to make nuclear shell model codes more transparent, easily accessible, and user friendly. Doesn't help that someone made a super performant code if it is not publicly available, difficult to obtain, and / or hard to use! Physicists are notorious when it comes to this, except maybe this one? 
Constant values and device memory
Changing the CPU calculation of the TBMEs to be GPU accelerated proved to be quite straight-forward actually. I had done most of the hard thinking when GPU accelerating the OBME calculation. The structure of the loops had to be slightly altered and I had to get rid of the Indices data structure inside the kernels. I solved this by chucking a bunch of constant device arrays in the global scope of hamiltonian_device.cpp. Global scope because I had a lot of trouble defining the constant device arrays in a separate file. The mentioned constant arrays contain correctly laid-out indices for the creation and annihilation operators in the shell model hamiltonian, as well as corresponding coupled $J$ and $M$ values. These values are pre-calculated in the ''generate_indices'' function and greatly speeds up the process of generating the TBMEs. All of these values are completely constant during the entire run time of the program and they are therefore aptly placed in the constant memory of the GPU (which apparently work a bit differently in AMD compared to Nvidia). The basis states are also constant during the entire run time of the program, but have to be placed in dynamic memory on the device because of the potentially huge size of the array. Well, they're not completely constant as they have to be calculated at the beginning of the program, but once they are calculated they do not change.
Discussion