From f2bb4a633d79a144887079ca73aa11da6ea50139 Mon Sep 17 00:00:00 2001
From: Niklaus Johner <nij2003@med.cornell.edu>
Date: Tue, 8 May 2012 11:23:35 -0400
Subject: [PATCH] Added function to calculate an RMSD matrix from a trajectory

---
 modules/mol/alg/pymod/trajectory_analysis.py | 38 +++++++++++++++++++-
 1 file changed, 37 insertions(+), 1 deletion(-)

diff --git a/modules/mol/alg/pymod/trajectory_analysis.py b/modules/mol/alg/pymod/trajectory_analysis.py
index 5e8c452b7..c23603205 100644
--- a/modules/mol/alg/pymod/trajectory_analysis.py
+++ b/modules/mol/alg/pymod/trajectory_analysis.py
@@ -4,7 +4,8 @@ Some functions for analyzing trajectories
 Author: Niklaus Johner
 """
 
-from ost import *
+import ost.mol.alg
+from ost import LogError
 import os
 
 def smooth(vec,n):
@@ -46,3 +47,38 @@ def smooth(vec,n):
   return vec2
 
 
+"""
+From here on the module needs numpy
+"""
+
+def RMSD_Matrix_From_Traj(t,sele,first=0,last=-1):
+  """
+  This function calculates a matrix M such that M[i,j] is the
+  RMSD of the EntityView sele between frames i and j of the trajectory t
+  aligned on sele.
+  Its inputs are:
+    t       : the trajectory (CoordGroupHandle)
+    sele    : the EntityView used for alignment and RMSD calculation
+    first=0 : the first frame of t to be used
+    last=-1 : the last frame of t to be used
+  Returns a numpy NxN matrix, where n is the number of frames.
+  """
+  try:
+    import numpy as npy
+    if last==-1:last=t.GetFrameCount()
+    n_frames=last-first
+    rmsd_matrix=npy.identity(n_frames)
+    for i in range(n_frames):
+      t=ost.mol.alg.SuperposeFrames(t,sele,begin=first,end=last,ref=i)
+      eh=t.GetEntity()
+      t.CopyFrame(i)
+      rmsd_matrix[i,:]=ost.mol.alg.AnalyzeRMSD(t,sele,sele)
+      if i==0:
+        last=last-first
+        first=0
+    return rmsd_matrix
+  except ImportError:
+    LogError("Function needs numpy, but I could not import it.")
+    raise
+
+
-- 
GitLab