Relax disp.spin lock offset+field figure

From relax wiki
Revision as of 00:13, 16 March 2014 by Tlinnet (talk | contribs)
Jump to navigation Jump to search

Reference to original figure

Comparing to Figure 1 and 10 the reference.

Palmer, A.G. & Massi, F. (2006). Characterization of the dynamics of biomacromolecules using rotating-frame spin relaxation NMR spectroscopy. Chem. Rev. 106, 1700-1719 | DOI: 

Script to produce figure

#############
# Made by Troels E. Linnet
# PhD student
# Copenhagen University
# SBiNLab, Structural Biology and NMR Laboratory.
# March 2014
#
# Trying to reproduce Figure 1 in:
#
# Palmer, A.G. & Massi, F. (2006). Characterization of the dynamics of biomacromolecules using rotating-frame spin relaxation NMR spectroscopy.
# Chem. Rev. 106, 1700-1719 http://dx.doi.org/10.1021/cr04042875
#
# This script is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#######################

import matplotlib.pyplot as plt
from matplotlib.patches import Circle, Ellipse, PathPatch, FancyArrowPatch
from mpl_toolkits.mplot3d import Axes3D, proj3d
import mpl_toolkits.mplot3d.art3d as art3d
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import pylab    

######## Define helper functions
class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)


def rotation_matrix(d):
    """
    Calculates a rotation matrix given a vector d. The direction of d
    corresponds to the rotation axis. The length of d corresponds to 
    the sin of the angle of rotation.

    Variant of: http://mail.scipy.org/pipermail/numpy-discussion/2009-March/040806.html
    """
    sin_angle = np.linalg.norm(d)

    if sin_angle == 0:
        return np.identity(3)

    d /= sin_angle

    eye = np.eye(3)
    ddt = np.outer(d, d)
    skew = np.array([[    0,  d[2],  -d[1]],
                  [-d[2],     0,  d[0]],
                  [d[1], -d[0],    0]], dtype=np.float64)

    M = ddt + np.sqrt(1 - sin_angle**2) * (eye - ddt) + sin_angle * skew
    return M

def pathpatch_2d_to_3d(pathpatch, z = 0, normal = 'z'):
    """
    Transforms a 2D Patch to a 3D patch using the given normal vector.

    The patch is projected into they XY plane, rotated about the origin
    and finally translated by z.
    """
    if type(normal) is str: #Translate strings to normal vectors
        index = "xyz".index(normal)
        normal = np.roll((1,0,0), index)

    normal /= np.linalg.norm(normal) #Make sure the vector is normalised

    path = pathpatch.get_path() #Get the path and the associated transform
    trans = pathpatch.get_patch_transform()

    path = trans.transform_path(path) #Apply the transform

    pathpatch.__class__ = art3d.PathPatch3D #Change the class
    pathpatch._code3d = path.codes #Copy the codes
    pathpatch._facecolor3d = pathpatch.get_facecolor #Get the face color    

    verts = path.vertices #Get the vertices in 2D

    d = np.cross(normal, (0, 0, 1)) #Obtain the rotation vector    
    M = rotation_matrix(d) #Get the rotation matrix

    pathpatch._segment3d = np.array([np.dot(M, (x, y, 0)) + (0, 0, z) for x, y in verts])

def pathpatch_translate(pathpatch, delta):
    """
    Translates the 3D pathpatch by the amount delta.
    """
    pathpatch._segment3d += delta

def update_position(e):
    print "From update position"
    #Transform co-ordinates to get new 2D projection
    tX, tY, _ = proj3d.proj_transform(dataX, dataY, dataZ, ax.get_proj())
    for i in range(len(dataX)):
        label = labels[i]
        label.xy = tX[i],tY[i]
        label.update_positions(fig.canvas.renderer)
    fig.canvas.draw()
    return

#########

fig = plt.figure()
ax=fig.gca(projection='3d')

# Draw a Center point
ax.scatter([0],[0],[0],c="k",s=50)

# Draw circle
circle = Circle((0, 0), 1)
circle.set_color("k")
circle.set_alpha(80)
#circle.set_fill(False)
ax.add_patch(circle)
#art3d.pathpatch_2d_to_3d(circle, z=0, zdir='y')
pathpatch_2d_to_3d(circle, z=0, normal = (1, 0, 1) )

#Input 3D Data for Text
#                  Sx        Sy       Sx         Szp             Sxp           w1        O1            we1        O_av          w_e
data = np.array([[1,0,0], [0,-1,0], [0,0,1], [0.975,0,1.17], [0.75,0,-0.75], [0.7,0,0], [0,0,0.3], [0.7,0,0.3], [0.0,0,0.84], [0.7,0,0.84]])
textlab = [r"S$_x$", r"S$_y$=S$_y$'", r"S$_z$", r"S$_z$'", r"S$_x$'", r"$\omega_1$", r"$\Omega_1$", r"$\omega_{e1}$", r"$\bar{\Omega}$", r"$\omega_{e}$"]

#Separate into X, Y, Z for greater clarity
dataX = data[:,0]
dataY = data[:,1]
dataZ = data[:,2]

#3D scatter plot
ax.scatter(dataX, dataY, dataZ, marker = 'o', c='k', s=1)

#Transform co-ordinates to get initial 2D projection
tX, tY, _ = proj3d.proj_transform(dataX, dataY, dataZ, ax.get_proj())

#Array of labels
labels = []

#Loop through data points to initially annotate scatter plot
#and populate labels array
for i in range(len(dataX)):
    #text='['+str(int(dataX[i]))+','+str(int(dataY[i]))+','+str(int(dataZ[i]))+']'
    text=textlab[i]
    label = ax.annotate(text,
            xycoords='data',
            xy = (tX[i], tY[i]), xytext = (-20, 20),
            textcoords = 'offset points', ha = 'right', va = 'top', fontsize=12,
            bbox = dict(boxstyle = 'round,pad=0.5', fc = 'grey', alpha = 0.8),
            arrowprops = dict(arrowstyle = '->', connectionstyle = 'arc3,rad=0'))
    labels.append(label)

# Make Sx
i = 0
Sx = Arrow3D([-1, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sx)

# Make Sy
i += 1
Sy = Arrow3D([0, dataX[i]],[1, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sy)

# Make Sz
i += 1
Sz = Arrow3D([0, dataX[i]],[0, dataY[i]], [-1, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(Sz)

# Make Szp
i += 1
Szp = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=3, arrowstyle="-|>", color="k")
ax.add_artist(Szp)

# Make Sxp
i += 1
Sxp = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=3, arrowstyle="-|>", color="k")
ax.add_artist(Sxp)

# Make w1
i += 1
w1 = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(w1)

# Make O1
i += 1
O1 = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(O1)

# Make we1
i += 1
we1 = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="b")
ax.add_artist(we1)

# Make O_av
i += 1
O_av = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-", color="k")
ax.add_artist(O_av)

# Make we
i += 1
we = Arrow3D([0, dataX[i]],[0, dataY[i]], [0, dataZ[i]], mutation_scale=20, lw=1, arrowstyle="-|>", color="r")
ax.add_artist(we)


ax.set_xlim3d(-1, 1)
ax.set_ylim3d(-1, 1)
ax.set_zlim3d(-1, 1)

fig.canvas.mpl_connect('button_release_event', update_position)
plt.show()