Recreate the projection of caffeine that you made in the first part. If you didn't make a script for this already, I suggest that you do so now. This file should create the hk-plane out to a given cut-off, calculate the molecular transform of caffeine on that plane, inverse transform to yield the projection, and display the projection with imshow.
Once you have recreated the projection, try to gradually decrease the cut-off. What happens to the atoms in the projection? What is the smallest cut-off that you can use and still distinguish individual (separated) atoms?
In the previous section, we mentioned that the resolution is the inverse of the maximum length of the scattering vector. Do you see what the relationship is between the maximum scattering vector and the cut-off variable? What is then the resolution you need to distinguish individual atoms?
Bragg's law [TODO: Add diffraction geometry diagram showing theta_max] One way to look at the concept of resolution is Bragg's law. In crystallography this is usually written as
In this equation θmax is half the maximum angle for which the scattered radiation is measured and dmin, the resolution, is interpreted as the smallest distance between resolvable lattice planes. Derived in slightly different way, the law says
Here, kmax is the maximum length of the scattering vector and kin is the length of the incident wave vector. From the left-hand side of the equations, you see that the resolution dmin is equivalent to 1/kmax. Do you see that the right-hand side of two equations are equivalent? (Remember the the definition of a wave vector). What does kmax correspond to in the planes that you created when calculating projections? In terms of the variables you used when creating projections, what is the resolution of the projections?
The convolution theorem Once you have accepted the relationship between diffraction and the Fourier transform of the electron density, it is often simpler to understand properties of diffraction in terms of the properties of the Fourier transform. For example, the electron density is a real function, and the Fourier transform of a real function is always hermitian. The absolute value of an hermitian function is centrosymmetric (symmetric around the origin, or even). Therefore, diffraction data, which represent the absolute square of the molecular transform, are always centrosymmetric. You can see this if you look at one of your planes in Fourier space.
We will now look at the resolution concept from the viewpoint of Fourier transforms. The whole story is expressed by one of the most important theorems in Fourier analysis, the convolution theorem:
The theorem says that if you take a function G, multiply by the function H and Fourier transform the result, you will end up with the Fourier transform of G (denoted g) convoluted with the Fourier transform of H (denoted h). The convolution theorem is exceptionally useful and one of the main reasons for the wide-spread use of Fourier transforms in computing today.
In the following section we will take a look at how the convolution theorem relates the resolution of an image to the distance that the diffraction data extends from the origin of Fourier space.
%%capture
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (8, 6)
# add path to the Modules
import sys, os
sys.path.append(os.path.join(os.getcwd(), '..', 'src'))
from Projection import *
# Read the pdb file for caffeine molecule.
caf = Molecule('caffeine.pdb')
# Distance between samples
step = 1/8.
cutoff = 10
H, K, L = TwoD_grid(step, cutoff)
F10 = moltrans(caf, H,K,L)
# Distance between samples
cutoff = 2
H, K, L = TwoD_grid(step, cutoff)
F2 = moltrans(caf, H,K,L)
rho10 = fft.fftshift(abs(fft.ifftn(F10, [2**10, 2**10])))
rho2 = fft.fftshift(abs(fft.ifftn(F2, [2**10, 2**10])))
W = squarewin2(F2.shape, F10.shape)
F10W = F10*W
rhoConv = fft.fftshift(abs(fft.ifftn(F10W, [2**10, 2**10])))
w = fft.fftshift(abs(fft.ifftn(W, [2**10, 2**10])))
Now plot F10, F2, F10W and W. Use subplot to display them all in a single figure. Set the colorbar range [-12,9] for F10, F2, F10W but for W to display the plots on the same color scale.
fig, ax = plt.subplots(2, 2)
ax[0,0].imshow(np.log(abs(F10)**2), interpolation='none',origin='lower', vmin=-12, vmax=9)
ax[0,0].set_title('F10')
ax[0,1].imshow(np.log(abs(F2)**2), interpolation='none',origin='lower', vmin=-12, vmax=9)
ax[0,1].set_title('F2')
tmp = abs(F10W)**2
tmp[tmp==0] = np.finfo(float).eps # replace 0 with epsilon to avoid problem with taking np.log.
ax[1,0].imshow(np.log(tmp), interpolation='none',origin='lower', vmin=-12, vmax=9)
ax[1,0].set_title('F10W')
ax[1,1].imshow(W, interpolation='none',origin='lower')
ax[1,1].set_title('W')
<matplotlib.text.Text at 0x11113f890>
Open a new figure and make four subplots containing rho10, rho2, rhoConv and w.
fig, ax = plt.subplots(2, 2)
ax[0,0].imshow(rho10, interpolation='none',origin='lower')
ax[0,1].imshow(rho2, interpolation='none',origin='lower')
ax[1,0].imshow(rhoConv, interpolation='none',origin='lower')
ax[1,1].imshow(w, interpolation='none',origin='lower')
<matplotlib.image.AxesImage at 0x103a46e10>
If you take a closer look at W, you will see that it has the same size as the plane F10 but that it takes value one within the area corresponding to F2 and the value zero outside. It follows that F2 = F10 · W and that the images F2 and F10W should be the identical but drawn at different scales. Check this by zooming in on F10W.
Now take a look at rhoConv. Applying what we said above and then the convolution theorem, we have
where w is the inverse Fourier transform of W. Hopefully, the images will show that rhoConv is identical to rho2.
The take-home message is the following: Since we don't measure the full, infinite, molecular transform (which we never can), we are in effect multiplying the infinite transform by some cut-off window (i.e. W). The effect of this in real space is to convolute the true density with the transform of the cut-off window. The effect of the convolution is to smear out the density, thus decreasing the resolution. In signal-processing terminology, w is called a low-pass lter. What is left after filtering is the low-frequency part of the density, corresponding to the inner part of the Fourier transform. W is called the transfer function of the filter.
Take another look at the Fourier transform of W. What you see is a two-dimensional sinc function,
The sinc function is also called the filtering function or interpolating function. It turns out that low-pass filtering and interpolation are related operations. If you look closely at the projection images, you can see the footprints of the sinc function in the form of electron density side-lobes around the atoms.
Plot the one-dimensional sinc function with Python.
x = np.linspace(-5, 5, 1001)
plt.plot(x, np.sinc(x))
[<matplotlib.lines.Line2D at 0x114f7de90>]
What is the width of the main lobe (central lobe) of the sinc-function? Try plot(x,sinc(2*x)). What is the width of the main lobe now? What is it for sinc(0.5x)? I guess this is pretty obvious from the equation. The interesting part comes here: The precise relationship between the sinc-function and the rectangle function W is
Or, in terms of the cutoff variable,
What is the relationship between the size of the cut-off in Fourier space, the width of the main lobe of the sinc function and what we call resolution?