If you are going to actually work through this tutorial, please work through examples/jupyter/quick start guide.ipynb first.
The cell below sets the notebook up for use with SpinEvolution. Run it each time after opening this notebook or after restarting the Jupyter kernel. If you are not going to run any other cells, there is no need to run this cell either.
#%matplotlib widget
%matplotlib inline
from plot import spinev, plot, setup_ipynb
setup_ipynb()
opts = '-v0 --verb --prev'
We attepmt to simulate ${\langle}I_{1z}+I_{2z}{\rangle}$ dephasing by the MELODRAMA recoupling pulse sequence as described in Sun et al. (1995), Internuclear distance measurements in solid state nuclear magnetic resonance: Dipolar recoupling via rotor synchronized spin locking, J.Chem.Phys. 102, 702-707.
We start with the main input file melo0, which we copy to melo:
%cp melo0 melo
Now open melo by double-clicking on it in Jupyter Lab File Browser. It will open up in a separate text editor window/tab. Drag this tab down to the middle of the screen and to the left to pace the Editor window side by side with the notebook. Then click on the File Browser tab to collapse the File Browser panel. We will keep this very convenient "split-screen" setup throughtout this whole tutorial, making edits to melo and using it as the main input file for the simulations.
The file used for each simulation can be always inspected by opening the collapsible panel right below the spinev() call. All the changes to this file (with respect to its version used in the pevious simulation) are shown automatically below the collapsible. These are the changes that you need to make (and save) in the Editor to follow through the Tutorial.
If something goes wrong, you can restart the kernel, run the setup cell to set everything up, and then proceed with the Tutorial starting from any place you want after copy-pasting the main input file used for that place from the drop-down panel into the Editor and saving it. Make sure that you do not copy the last line (options_cmd_line).
We are ready for our first simulation now:
spinev('melo',opts)
The results of the similation look very different from Fig. 2 of the paper.
How can that be?
Let's first check out the pulse sequence:
spinev('melo -pseq',opts)
The pulse sequence looks good: we see that it is indeed $R\overline{R}\overline{R}R$, where $R = XYYX$
Let's look at the magnetization path in more detail, observing it after every RF step during the whole evolution time. This is done by using the -oes option and making the pulse sequence zero-dimensional with 10 repeats.
spinev('melo -oes',opts)
And now let's make sure that this suspiciouly wild-looking curve indeed contains all the points from the curve we simulated earlier: we need to be sure that we get identical results using these two different approaches. The first approach was simple 1D pulse sequence; the second – 0D pulse sequence supplemented with the -oes option.
In order to do that, we revert the pulse sequence back to 1D and run the simulation again, but now saving the output under a different name so that it does not overwrite the output of our last simulation (which was set up as 0D+oes), and using this previuos output as a reference on the plot.
This can be done with the help of the options -n<outname>) and -ref<refname>, respectively. We also use --lines and --markers plot options to plot the output as small filled circles, and the reference – as the usual solid line.
To show the changes in the main input file, we use plot option --ref<refname> because the output is saved under a different name.
spinev('melo -ntmp --lines=" ",- --markers=." " -refmelo --ref melo',opts)
Indeed, everything is coorect.
Let's try to understand how exactly this "wild" curve comes about. It is a result of several things going on at the same time. There is flipping of the magnetization with every RF pulse. There is dipolar evolution. And there is powder averaging.
Quite obviously, the curve is jagged due to the RF pulses. It would be nice if we could simplify their effects. Recall that each pulse is $\tau_R/4$ long and its power is N times $1/\tau_R$. Which means that its flip angle is N/4 full turns. So if N is a multiple of 4, this jaggedness should dissapear!
Let's try it.
spinev('melo -oes --ref tmp',opts)
Yes!
But now we have another problem. This dephasing curve is obviously incorrect. And this kind of "incorrectness" is usually due to insufficient powder averaging. We need to either increase the number of $(\alpha,\beta)$ angles in the set, or the number of steps in the averaging over $\gamma$, or both. Let's try gamma at first because 16 angles, althouth fast to calculate, may be a very poor choice due to the pulse sequence symmetry.
spinev('melo -oes',opts)
Looks good for now.
Let's now repeat the 1D vs. 0D+oes comparison for N=4:
spinev('melo -ntmp --lines=" ",- --markers=." " -refmelo --ref melo',opts)
Looks right. But our main problem is still not resolved: N=4 and N=5 curves are dramatically different. Let's look at them on the same plot:
spinev('melo -oes --ref tmp',opts)
The N=4 curve decays much faster than the N=5 one. How is it possible? The only difference in these two simulations is the 25% difference in the RF power. According to Sun et al., MELODRAMA design implies that the dephasing curves for these two cases will be very similar.
To understand this, let's look at the spin dynamics in maximum detail. We will start with removing the powder averaging:
spinev('melo -oes',opts)
The dephasing of ${\langle}I_{1z}+I_{2z}{\rangle}$ occurs because each crystallite has a different frequency of this oscillation. And in our case, the frequencies for N=4 and N=5 are indeed very simular, as should be expected for MELODRAMA. The problem then must be in how these frequencies are modulated, i.e. in this jagged pattern produced by the RF.
Let's look then at this modulation in pure form, i.e. without the dipolar interaction.
Option -scheck13 simply suppresses the warning that appears when we attempt to simulate evolution under isotropic Hamiltonian while using MAS.
spinev('melo -oes --markers=. -scheck13',opts)
To see the full detail, zoom in on the first two RF cycles:
spinev('melo -oes --markers=. -scheck13',opts)
This looks correct: Each RF pulse in the N=4 case is a $2\pi$ rotation, while the N=5 pulses are $\pi/2$ rotations whith phases $XYYX$ or $\overline{XYYX}$.
Let's restore the dipolar interaction now:
spinev('melo -oes --markers=.',opts)
Now, this is the moment of truth.
We see that the RF behaves exactly as specified by the MELODRAMA pulse sequence. But the evolution due to the diploar Hamiltonian restored by this pulse sequence is dramatically different for N=4 and N=5 cases even on the time scale of pulse sequence cycle. Which means that the average Hamiltonians produced by these two cases are very different. Which means that what the paper states about this pulse sequence is wrong.
Now, depending on how your intuition works, there can be many ways to guess was is right then. Make your own story here.
My story was as follows. I was looking at the picture above and marveling at how messed up this $XYYX$ phase cycle is. It is not even cyclic for N=5, i.e. it does not return the magnetization to its original state after one cycle. There certainly must be a better way of doing this kind of recoupling. And then I remembered that Fig. 1 caption of the paper mentioned a modification of this pulse sequence that "partially compensates for rf field inhomogeneity", and that this modification was $X\overline{X}Y\overline{Y}\overline{Y}Y\overline{X}X$. Well, but these $X\overline{X}$ and $Y\overline{Y}$ combinations are exactly what we need to make this phase cycle right! Let's implement it!
Since there are now 32 steps in the RF cycle, it will be quicker and more reliable if we compute it rather than create it in a text editor. So we will use the [text]N macro instead of melo.pp to define the pulse sequence. The macro simply expands to N space-separated copies of text. We will also use the modulo function to compute $\overline{R}$ from $R$, and the ; operator to concatenate $R$ and $\overline{R}$ column vectors into the full supercycle:
spinev('melo -oes --markers=.',opts)
Now the N=4 and N=5 cases behave very similarly: our guess was correct!
It remains now to sample the data properly, i.e. once every rotor period or something like that. If we simulate the pulse sequence as 1D, we are limited to sampling it once per full RF cycle, i.e. once every $4\tau_R$. So we'll stick with our 0D+oes approach.
The size of the output matrix, pdata_re, cannot be changed and we cannot somehow "skip" points from it. So we will leave the standard output as is, and will simply supplement it with our own output, which we will create with the plot() function available in the Variables section.
We use getind(data_label>=0) to get all element indices of the data_label vector, which consists of the observation times for each point. This works because the condition data_label>=0 is True for all points. The next line (see below) filters out every 8th point with the help of the modulo function.
Note that there are two plots on the output: one poroduced by our plot() call, and one that visualizes the standard simulation output.
spinev('melo -oes --markers=.',opts)
Now restore the powder averaging and the full evolution length:
spinev('melo -oes',opts)
And let's see if we get a smoother curve if we sample twice per rotor cycle:
spinev('melo -oes',opts)
Looks perfect now.
So let's get the curves for other N:
spinev('melo -oes',opts)
spinev('melo -oes',opts)
Now let's go back to the basic MELODRAMA sequence and compute the dephasing for odd and even N separately:
spinev('melo -oes',opts)
spinev('melo -oes',opts)
We can see that the basic (4-step) MELODRAMA sequence has three cases for the recoupled Hamiltonian it generates: N=1, odd N, and even N.
The cyclic (8-step) version of MELODRAMA recouples more uniformly for different N, except for N=1. The strongest recoupling is observed for the N=3 case.
The analytic (first-order AHT) estimate of the recoupled Hamiltonian is off by at least a factor of 2 for the N=1 case with both versions of the sequence.
It appears that the authors did not realize that the basic version of their pulse sequence has very different average Hamiltonians for odd and even N. The 8-step version is also somewhat sensitive to the N value, but to a much lesser degree. Thus, experimantally, it may appear that it just compensates for RF inhomogeneity, while in reality it produces a very different recoupling.