In [1]:

    

# In this tutorial we'll go over how to simulate
# x-ray scattering on a polar detector. This is
# a fast way to test out what correlations should
# look like.

from odin import xray
from mdtraj import trajectory


    

In [2]:

    

# load in a PDB file
t = trajectory.load('3LYZ.pdb') # your fav pdb


    

In [25]:

    

n_shots = 5                         # total number of shots to do
n_molecules = 5                     # the number of molecules to include per shot
q_values = np.arange(1.2,1.6,0.02)  # the |q| values of the rings to sim
n_phi = 360                         # number of pts around the rings

rings = xray.Rings.simulate(t, n_molecules, q_values, n_phi, n_shots)


    

    




15:25:53 - Running 5 molecules from snapshot 0...
15:25:55 - Finished polar shot 1/5 on device 0
15:25:55 - Running 5 molecules from snapshot 0...
15:25:57 - Finished polar shot 2/5 on device 0
15:25:57 - Running 5 molecules from snapshot 0...
15:25:58 - Finished polar shot 3/5 on device 0
15:25:58 - Running 5 molecules from snapshot 0...
15:26:00 - Finished polar shot 4/5 on device 0
15:26:00 - Running 5 molecules from snapshot 0...
15:26:02 - Finished polar shot 5/5 on device 0

In [26]:

    

# check the WAXS
figure()

waxs = rings.intensity_profile() # first col is |q|, second is intensity

plot(waxs[:,0], waxs[:,1], lw=2)
xlabel(r'$q / \AA^{-1}$')
ylabel('Intensity')
show()


    

In [27]:

    

# looks like we have a peak at |q| = 1.46, so let's correlate that
# I'll make a plot of each ring's correlation function independently

corr = rings.correlate_intra(1.46, 1.46) # args are q_1, q_2

figure()

for i in range(rings.num_shots):
    plot(rings.phi_values, corr[i,:], lw=2)

xlabel(r'$\phi$ / Radians')
ylabel(r'$C(\phi)$')
show()


    

In [28]:

    

# let's check out the average of each of those shots
# also include a visualization of the std

figure()
plot(rings.phi_values, corr.mean(0), lw=2)
plot(rings.phi_values, corr.mean(0) + corr.std(0), lw=0.5, color='k')
plot(rings.phi_values, corr.mean(0) - corr.std(0), lw=0.5, color='k')
xlabel(r'$\phi$ / Radians')
ylabel(r'$C(\phi)$')
show()


    

In [ ]: