hdiff output

r27332/hbond_matrix.py 2015-11-17 23:33:58.949295311 +0000 r27331/hbond_matrix.py 2015-11-17 23:33:59.157298103 +0000
  1: #!/usr/bin/python  1: #!/usr/bin/python
  2:   2: 
  3: def make_hbond_input(hbond_input_filename, hbond_output_filename,  3: def make_hbond_input(hbond_input_filename, hbond_output_filename,
  4:                      donor_acceptor_filename, restart_filename,  4:                      donor_acceptor_filename, restart_filename, 
  5:                      distance_cutoff, angle_cutoff):  5:                      distance_cutoff, angle_cutoff):
  6:     """ 
  7:     Make the hbond input file with specified cutoffs and input/output files. 
  8:     """ 
  9:     with open(hbond_input_filename, 'w') as hbond_input:  6:     with open(hbond_input_filename, 'w') as hbond_input:
 10:         # trajin command to read restart file  7:         # trajin command        
 11:         hbond_input.write(' '.join(['trajin', restart_filename, '\n']))  8:         hbond_input.write(' '.join(['trajin', restart_filename, '\n']))
 12:         # contents of donor/acceptor file to specify hbonding groups  9:         # contents of donor/acceptor file
 13:         with open(donor_acceptor_filename, 'r') as donor_acceptor_file: 10:         with open(donor_acceptor_filename, 'r') as donor_acceptor_file:
 14:             for line in donor_acceptor_file: 11:             for line in donor_acceptor_file:
 15:                 hbond_input.write(line) 12:                 hbond_input.write(line)
 16:         # trajout command to write output file 13:         # trajout command
 17:         hbond_input.write(' '.join(['hbond', 'out', hbond_output_filename, 14:         hbond_input.write(' '.join(['hbond', 'out', hbond_output_filename, 
 18:                                     'distance', str(distance_cutoff), 15:                                     'distance', str(distance_cutoff), 
 19:                                     'angle', str(angle_cutoff), '\n'])) 16:                                     'angle', str(angle_cutoff), '\n']))
 20:  17: 
 21: def run_ptraj(hbond_input_filename, topology_filename): 18: def run_ptraj(hbond_input_filename, topology_filename):
 22:     """ 
 23:     Run ptraj as a subprocess, using the hbond input file as STDIN and piping 
 24:     STDOUT to /dev/null. 
 25:     """ 
 26:     import os 19:     import os
 27:     import os.path 20:     import os.path
 28:     import subprocess 21:     import subprocess
 29:     with open(hbond_input_filename, 'r') as hbond_input: 22:     with open(hbond_input_filename, 'r') as hbond_input:
 30:         with open(os.devnull, 'w') as devnull: 23:         with open(os.devnull, 'w') as devnull:
 31:             # Runs ptraj 
 32:             ptraj = subprocess.Popen([os.path.join(os.environ['AMBERHOME'], 'bin', 'ptraj'), topology_filename], stdin=hbond_input, stdout=devnull) 24:             ptraj = subprocess.Popen([os.path.join(os.environ['AMBERHOME'], 'bin', 'ptraj'), topology_filename], stdin=hbond_input, stdout=devnull)
 33:             ptraj.wait() 25:             ptraj.wait()
 34:  26: 
 35: def get_residues(residues_filename, topology_filename): 27: def get_residues(residues_filename, topology_filename):
 36:     """ 
 37:     Get the residues from the residue file, if specified. If the file is 
 38:     not specified, use the list of all residues from the topology file. 
 39:     """ 
 40:     residues = [] 28:     residues = []
 41:     # If residue file is specified, use those residues. 29:     if residues_filename: 
 42:     if residues_filename: 
 43:         with open(residues_filename, 'r') as residue_input: 30:         with open(residues_filename, 'r') as residue_input:
 44:             for line in residue_input: 31:             for line in residue_input:
 45:                 residues.append(int(line.strip())) 32:                 residues.append(int(line.strip()))
 46:     # Otherwise, use all residues from the topology file. 
 47:     else: 33:     else:
 48:         with open(topology_filename, 'r') as prmtop: 34:         with open(topology_filename, 'r') as prmtop:
 49:             target_line_num = None 35:             target_line_num = None
 50:             for i, line in enumerate(prmtop): 36:             for i, line in enumerate(prmtop):
 51:                 if 'POINTERS' in line: 37:                 if 'POINTERS' in line:
 52:                     target_line_num = i + 3 38:                     target_line_num = i + 3
 53:                 if i == target_line_num: 39:                 if i == target_line_num:
 54:                     num_residues = int(line.split()[1]) 40:                     num_residues = int(line.split()[1])
 55:                     residues = range(1, num_residues + 1) 41:                     residues = range(1, num_residues + 1)
 56:     return residues 42:     return residues
 57:  43: 
 58: def create_hbond_dict(hbond_output_filename): 44: def create_hbond_dict(hbond_output_filename): 
 59:     """ 45:     # hbond_dict keys are (donor_res, 'bb'|'sc', acc_res, 'bb'|'sc')
 60:     Creates a dictionary containing the hydrogen bonds found by ptraj. This has 46:     # hbond_dict values are the number of bonds 
 61:     keys which are tuples of the form: 
 62:  
 63:     (donor_res, 'bb'|'sc', acc_res, 'bb'|'sc') 
 64:      
 65:     The values are the number of bonds. 
 66:      
 67:     e.g. 
 68:         hbond_dict[(12, 'bb', 15, 'sc')] = 2 indicates that residue 12 is a  
 69:         backbone donor to the sidechain of residue 15. 
 70:     """ 
 71:     hbond_dict = {} 47:     hbond_dict = {}
 72:     with open(hbond_output_filename, 'r') as hbond_output: 48:     with open(hbond_output_filename, 'r') as hbond_output:
 73:         for line in hbond_output: 49:         for line in hbond_output:
 74:             words = line.split() 50:             words = line.split()
 75:             try: 51:             try:
 76:                 if words[0] == '|': 52:                 if words[0] == '|':
 77:                     donor_res, donor_atom = map(str.strip, words[2].strip(':').split('@')) 53:                     donor_res, donor_atom = map(str.strip, words[2].strip(':').split('@'))
 78:                     acceptor_res, acceptor_atom = map(str.strip, words[7].strip(':').split('@')) 54:                     acceptor_res, acceptor_atom = map(str.strip, words[7].strip(':').split('@'))
 79:                     # Make sure that the residue numbers are integers. 55:                     # Make sure that the residue numbers are integers.
 80:                     donor_res, acceptor_res = map(int, (donor_res, acceptor_res)) 56:                     donor_res, acceptor_res = map(int, (donor_res, acceptor_res))
 88:                         acceptor_type = 'sc' 64:                         acceptor_type = 'sc'
 89:                     try: 65:                     try:
 90:                         hbond_dict[(donor_res, donor_type, acceptor_res, acceptor_type)] += 1 66:                         hbond_dict[(donor_res, donor_type, acceptor_res, acceptor_type)] += 1
 91:                     except KeyError: 67:                     except KeyError:
 92:                         hbond_dict[(donor_res, donor_type, acceptor_res, acceptor_type)] = 1 68:                         hbond_dict[(donor_res, donor_type, acceptor_res, acceptor_type)] = 1
 93:             except IndexError: 69:             except IndexError:
 94:                 continue 70:                 continue
 95:     return hbond_dict 71:     return hbond_dict
 96:  72: 
 97: def print_matrix(hbond_dict, residues): 73: def print_matrix(hbond_dict, residues):
 98:     """ 
 99:     Prints the hbond dictionary as a matrix with two rows per donor residue 
100:     (one backbone, one sidechain) and two columns per acceptor residue. 
101:     """ 
102:     for res1 in residues: 74:     for res1 in residues:
103:         for res1_type in ['bb', 'sc']: 75:         for res1_type in ['bb', 'sc']:
104:             for res2 in residues: 76:             for res2 in residues:
105:                 for res2_type in ['bb', 'sc']: 77:                 for res2_type in ['bb', 'sc']:
106:                     try: 78:                     try:
107:                         print hbond_dict[(res1, res1_type, res2, res2_type)], 79:                         print hbond_dict[(res1, res1_type, res2, res2_type)],
108:                     except KeyError: 80:                     except KeyError:
109:                         print "0", 81:                         print "0",
110:             print "" 82:             print ""
111:  83: 
115:     # Assign the command-line arguments to specify file names 87:     # Assign the command-line arguments to specify file names
116:     topology_filename, restart_filename, donor_acceptor_filename = sys.argv[1:4] 88:     topology_filename, restart_filename, donor_acceptor_filename = sys.argv[1:4]
117:     # If specified, get a list of residues we're interested in 89:     # If specified, get a list of residues we're interested in
118:     if len(sys.argv) >= 5: 90:     if len(sys.argv) >= 5:
119:         residues_filename = sys.argv[4] 91:         residues_filename = sys.argv[4]
120:     else: 92:     else:
121:         residues_filename = None 93:         residues_filename = None
122:     # If specified, set a custom distance and angle cutoff 94:     # If specified, set a custom distance and angle cutoff
123:     if len(sys.argv) == 7: 95:     if len(sys.argv) == 7:
124:         distance_cutoff, angle_cutoff = map(float, sys.argv[5:7]) 96:         distance_cutoff, angle_cutoff = map(float, sys.argv[5:7])
125:     # Otherwise use defaults of 3.0A and 120 degrees. 
126:     else: 97:     else:
127:         distance_cutoff, angle_cutoff = 3.0, 120.0 98:         distance_cutoff, angle_cutoff = 3.0, 120.0
128:  99: 
129:     # Set default names for the hbond files100:     # Set default names for the hbond files
130:     hbond_input_filename, hbond_output_filename = 'hbond.in', 'hbond.out'101:     hbond_input_filename, hbond_output_filename = 'hbond.in', 'hbond.out'
131: 102: 
132:     # Make the input file103:     # Make the input file
133:     make_hbond_input(hbond_input_filename, hbond_output_filename,104:     make_hbond_input(hbond_input_filename, hbond_output_filename,
134:                      donor_acceptor_filename, restart_filename,105:                      donor_acceptor_filename, restart_filename,
135:                      distance_cutoff, angle_cutoff)106:                      distance_cutoff, angle_cutoff)
136: 107: 
137:     # Run ptraj108:     # Run ptraj    
138:     run_ptraj(hbond_input_filename, topology_filename)109:     run_ptraj(hbond_input_filename, topology_filename)
139: 110:    
140:     # Get the list of important residues111:     # Get the list of important residues
141:     residues = get_residues(residues_filename, topology_filename)112:     residues = get_residues(residues_filename, topology_filename)
142: 113: 
143:     # Create the dictionary corresponding to the hbond matrix114:     # Create the dictionary corresponding to the hbond matrix
144:     hbond_dict = create_hbond_dict(hbond_output_filename)115:     hbond_dict = create_hbond_dict(hbond_output_filename)
145: 116: 
146:     # Print the matrix117:     # Print the matrix
147:     print_matrix(hbond_dict, residues)118:     print_matrix(hbond_dict, residues)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0