Solutions to the exercises of part 4

Functions of a matrix

import Numeric; N = Numeric
import LinearAlgebra; LA = LinearAlgebra

# Calculate an arbitrary function of a matrix
def matrixFunction(function, matrix):
    eigenvalues, eigenvectors = LA.eigenvectors(matrix)
    eigenvalues = function(eigenvalues)
    return N.dot(eigenvalues*N.transpose(eigenvectors), eigenvectors)

# Calculate an exponential by Taylor expansion
def matrixExponential(matrix):
    sum = N.identity(matrix.shape[0])
    product = sum
    for i in range(1, 50):
	product = N.dot(product, matrix)/i
	sum = sum + product
    return sum

# Create a symmetric matrix.
m = N.reshape(N.arrayrange(9.), (3, 3))**2/10.
m = (m + N.transpose(m))/2.

# Compare the results
print matrixExponential(m)
print matrixFunction(N.exp, m)

Visualization of differences in conformations

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.Geometry import Vector
import string
from Scientific.Visualization import VRML

generic_format = FortranFormat('A6')
atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

# Read the PDB file and make a list of all C-alpha positions
def readCAlphaPositions(filename):
    c_alpha_positions = []
    for line in TextFile(filename):
	record_type = FortranLine(line, generic_format)[0]
	if record_type == 'ATOM  ' or record_type == 'HETATM':
	    data = FortranLine(line, atom_format)
	    atom_name = string.strip(data[2])
	    position = Vector(data[8:11])
	    if atom_name == 'CA':
		c_alpha_positions.append(position)
    return c_alpha_positions

# Read the two conformations
conf1 = readCAlphaPositions('conf1.pdb')
conf2 = readCAlphaPositions('conf2.pdb')

# Create the VRML scene and the material for the arrows
scene = VRML.Scene([])
arrow_material = VRML.DiffuseMaterial('black')
arrow_radius = 0.2

# Create the spheres and put them into the scene
for i in range(len(conf1)):
    scene.addObject(VRML.Arrow(conf1[i], conf2[i], arrow_radius,
			       material=arrow_material))

# View the scene
scene.view()