folder name change and README minor format

interior_residue/add_zeros_to_fnm.py

@@ -0,0 +1,71 @@
+import numpy as np
+import sys
+import os
+#from math import exp
+#from Bio.PDB import *
+
+
+'''
+For each FNM file -- add 0-vectors to so that 10 modes are avail
+
+USAGE:
+python add_zeros_to_fnm.py 
+stdin/stdout
+
+'''
+
+
+
+def check_if_list_has_only_floats(list_to_check):
+	status = True
+	list_of_float_like_chars = ("-", ".", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9")
+	for element in list_to_check:
+		for charact in element:
+			if charact not in list_of_float_like_chars:
+				return False
+	return status
+
+
+num_modes = 0
+avail_fnm_file_lines = list()
+fnm_file_to_rd = sys.stdin
+mode___2___num_vectors = {}
+number_of_vectors_in_mode = 0  ## re-setting the counter
+for line in fnm_file_to_rd:
+	avail_fnm_file_lines.append(line)
+	if line[0] is not "#":  ## ie, if we're not at file header
+		ln_elmnts = list()
+		ln_elmnts = line.split()
+		if "BEGIN MODE" in line:
+			num_modes += 1
+			mode_indx = int(ln_elmnts[2])
+			number_of_vectors_in_mode = 0  ## re-setting the counter
+		if (len(ln_elmnts) == 3)  and  (check_if_list_has_only_floats(ln_elmnts)):
+			number_of_vectors_in_mode += 1
+		if line[0:3] == "END":
+			mode___2___num_vectors[mode_indx] = number_of_vectors_in_mode
+fnm_file_to_rd.close()
+	#for mode in mode___2___num_vectors:
+	#	print fnm_file + "    mode: " + str(mode) + "     num_modes: " + str(num_modes)
+
+out_file_name_to_wrt = sys.stdout
+for avail_line in avail_fnm_file_lines:
+	out_file_name_to_wrt.write(avail_line)
+
+
+num_remaining_modes_to_print = 10 - num_modes
+starting_mode = 17 - num_remaining_modes_to_print
+mode_to_print = starting_mode
+while mode_to_print <= 16:
+	out_file_name_to_wrt.write("BEGIN MODE " + str(mode_to_print) + "\n")
+	vect = 0
+	while vect < number_of_vectors_in_mode:
+		out_file_name_to_wrt.write("0.0000000000 0.0000000000 0.0000000000\n")
+		vect += 1
+	out_file_name_to_wrt.write("END\n\n")
+	mode_to_print += 1
+out_file_name_to_wrt.close()
+
+
+#print "\n\n --- run complete -- \n\n"
+

interior_residue/calpha_modes.py

@@ -0,0 +1,106 @@
+#!python
+
+# This example shows how to calculate approximate low-frequency
+# modes for big proteins. For a description of the techniques,
+# see
+#
+#    K. Hinsen
+#    Analysis of domain motions by approximate normal mode calculations
+#    Proteins 33, 417 (1998)
+#
+# and
+#
+#    K. Hinsen, A.J. Petrescu, S. Dellerue, M.C. Bellissent-Funel, G.R. Kneller
+#    Harmonicity in slow protein dynamics
+#    Chem. Phys. 261, 25 (2000)
+#
+##
+#
+# Adaptation by Simon Mitternacht 2010-2011.
+#
+
+
+from MMTK import *
+from MMTK.Proteins import Protein
+from MMTK.Proteins import PeptideChain
+from MMTK.ForceFields import CalphaForceField
+from MMTK.FourierBasis import FourierBasis, estimateCutoff
+from MMTK.NormalModes import NormalModes, VibrationalModes, SparseMatrixSubspaceNormalModes
+from numpy import *
+import sys
+import re
+
+
+# top_nm.pl only reads the first 20 modes before exiting. This results in a broken pipe error.
+# Prevent calpha_modes.py from crashing by ignoring SIGPIPE signal.
+import signal
+signal.signal(signal.SIGPIPE, signal.SIG_DFL)
+
+
+
+def modes_in_sequential_order(protein,modes):
+    modes_inseq = []
+    for mode in modes: #[modes.rawMode(i) for i in xrange(len(modes))]:
+        tmp_mode = []
+        for chain in protein:
+            for residue in chain.residues():
+                tmp = mode[residue.peptide.C_alpha]
+                norm = 1#math.sqrt(residue.mass())
+                tmp_vec = [tmp.x()/norm, tmp.y()/norm, tmp.z()/norm]
+                tmp_mode.append(tmp_vec)
+        modes_inseq.append(tmp_mode)
+    return modes_inseq
+
+
+
+if (len(sys.argv) != 2 and len(sys.argv) != 3):
+    print "Usage: %s pdb [v|f]" % sys.argv[0]
+    print "Outputs the 10 % lowest frequency normal modes,to stdout. "
+    print "Second optional argument specifies if VibrationalModes or FourierBasis should be used."
+    exit()
+
+
+# Check CL arguments    
+pdb_file = sys.argv[1]
+type = "auto"
+if (len(sys.argv) > 2):
+    arg = sys.argv[2]
+    if (arg != 'v' and arg != 'f'):
+        sys.stderr.write("Argument \'%s\' invalid" % arg)
+        exit()
+    type = arg
+
+# Construct system
+universe = InfiniteUniverse(CalphaForceField(2.5))
+universe.protein = Protein(pdb_file, model='calpha')
+
+# Find a reasonable basis set size and cutoff
+cutoff = None
+if (type == "auto") or (type == "f"):
+    nbasis = max(10, universe.numberOfAtoms()/10)
+    cutoff, nbasis = estimateCutoff(universe, nbasis)
+
+# calculate modes for first system
+if (cutoff is None):
+    # Do full normal mode calculation
+    print "# VibrationalModes"
+    modes = VibrationalModes(universe)
+    print "# %d" % len(modes)
+else:
+    # Do subspace mode calculation with Fourier basis for large proteins
+    print "# FourierBasis: cutoff = %d, nbasis = %d" % (cutoff, nbasis)
+    subspace = FourierBasis(universe, cutoff)
+    subspace.may_modify = 1
+    modes = SparseMatrixSubspaceNormalModes(universe, subspace)
+
+modes_seq = modes_in_sequential_order(universe.protein, modes)
+N = min(len(modes_seq),3*universe.numberOfAtoms()/10+6)
+
+for i in range(6,N):
+    print "BEGIN MODE %d" % ((i+1))
+    for d in modes_seq[i]:
+        print "%.10f %.10f %.10f" % (d[0], d[1], d[2])
+    print "END"
+    print ""
+
+