diff --git a/modules/mol/alg/pymod/trajectory_analysis.py b/modules/mol/alg/pymod/trajectory_analysis.py index c236032051edd9eee2069f02c26eb3f5b71f5532..f71e42740657d4c7a7d750403342ea86f8e317b8 100644 --- a/modules/mol/alg/pymod/trajectory_analysis.py +++ b/modules/mol/alg/pymod/trajectory_analysis.py @@ -82,3 +82,68 @@ def RMSD_Matrix_From_Traj(t,sele,first=0,last=-1): raise +def PairwiseDistancesFromTraj(t,sele,first=0,last=-1,seq_sep=1): + """ + This function calculates the distances between any pair of atoms in the + EntityView sele with sequence separation larger than seq_sep from a trajectory t. + It return a matrix containing one line for each atom pair and N columns, where + N is the length of the trajectory. + Its inputs are: + t : the trajectory (CoordGroupHandle) + sele : the EntityView used to determine the atom pairs + first=0 : the first frame of t to be used + last=-1 : the last frame of t to be used + seq_sep=1 : The minimal sequence separation between + Returns a numpy NpairsxNframes matrix. + """ + try: + import numpy as npy + if last==-1:last=t.GetFrameCount() + n_frames=last-first + n_var=0 + for i,a1 in enumerate(sele.atoms): + for j,a2 in enumerate(sele.atoms): + if not j-i<seq_sep:n_var+=1 + #n_var=sele.GetAtomCount() + #n_var=(n_var-1)*(n_var)/2. + dist_matrix=npy.zeros(n_frames*n_var) + dist_matrix=dist_matrix.reshape(n_var,n_frames) + k=0 + for i,a1 in enumerate(sele.atoms): + for j,a2 in enumerate(sele.atoms): + if j-i<seq_sep:continue + dist_matrix[k]=ost.mol.alg.AnalyzeDistanceBetwAtoms(t,a1.GetHandle(),a2.GetHandle())[first:last] + k+=1 + return dist_matrix + except ImportError: + LogError("Function needs numpy, but I could not import it.") + raise + +def DistanceMatrixFromPairwiseDistances(distances,p=2): + """ + This function calculates an distance matrix M(NxN) from + the pairwise distances matrix D(MxN), where N is the number + of frames in the trajectory and M the number of atom pairs. + M[i,j] is the distance between frame i and frame j + calculated as a p-norm of the differences in distances + from the two frames (distance-RMSD for p=2). + Inputs: + distances : a pairwise distance matrix as obtained from PairwiseDistancesFromTraj() + Returns a numpy NxN matrix, where N is the number of frames. + """ + try: + import numpy as npy + n1=distances.shape[0] + n2=distances.shape[1] + dist_mat=npy.identity(n2) + for i in range(n2): + for j in range(n2): + if j<=i:continue + d=(((abs(distances[:,i]-distances[:,j])**p).sum())/float(n1))**(1./p) + dist_mat[i,j]=d + dist_mat[j,i]=d + return dist_mat + except ImportError: + LogError("Function needs numpy, but I could not import it.") + raise +