Coverage for model_workflow/utils/vmd_spells.py: 79%
292 statements
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-23 10:54 +0000
« prev ^ index » next coverage.py v7.9.2, created at 2025-07-23 10:54 +0000
1# Functions powered by VMD
2# Humphrey, W., Dalke, A. and Schulten, K., "VMD - Visual Molecular Dynamics", J. Molec. Graphics, 1996, vol. 14, pp. 33-38.
3# http://www.ks.uiuc.edu/Research/vmd/
5import os
6from os.path import exists
8from subprocess import run, PIPE, STDOUT
10from model_workflow.utils.file import File
11from model_workflow.utils.type_hints import *
12from model_workflow.utils.auxiliar import warn
14# Set characters to be escaped since they have a meaning in TCL
15TCL_RESERVED_CHARACTERS = ['"','[',']']
17# Given a VMD atom selection string, escape TCL meaningful characters and return the escaped string
18def escape_tcl_selection (selection : str) -> str:
19 escaped_selection = selection
20 for character in TCL_RESERVED_CHARACTERS:
21 escaped_selection = escaped_selection.replace(character, '\\' + character)
22 return escaped_selection
24# Set the script filename with all commands to be passed to vmd
25commands_filename = '.commands.vmd'
27# List all the vmd supported trajectory formats
28vmd_supported_structure_formats = {'pdb', 'prmtop', 'psf', 'parm', 'gro'} # DANI: Esto lo he hecho rápido, hay muchas más
29vmd_supported_trajectory_formats = {'mdcrd', 'crd', 'dcd', 'xtc', 'trr', 'nc', 'netcdf', 'cdf', 'pdb', 'rst7'}
31# Set a vmd format translator
32vmd_format = {
33 'pdb': 'pdb',
34 'prmtop': 'prmtop',
35 'psf': 'psf',
36 'gro': 'gro',
37 'crd': 'crd',
38 'mdcrd': 'crd',
39 'dcd': 'dcd',
40 'xtc': 'xtc',
41 'trr': 'trr',
42 'netcdf': 'netcdf',
43 'nc': 'netcdf',
44 'cdf': 'netcdf',
45 'rst7': 'rst7'
46}
48# Given a vmd supported topology with no coordinates and a single frame file, generate a pdb file
49def vmd_to_pdb (
50 input_structure_filename : str,
51 input_trajectory_filename : str,
52 output_structure_filename : str):
54 # Get the trajectory format according to the vmd dictionary
55 input_trajectory_file = File(input_trajectory_filename)
56 trajectory_format = input_trajectory_file.format
57 vmd_trajectory_format = vmd_format[trajectory_format]
59 # Prepare a script for VMD to run. This is Tcl language
60 with open(commands_filename, "w") as file:
61 # Load only the first frame of the trajectory
62 file.write(f'animate read {vmd_trajectory_format} {input_trajectory_filename} end 1\n')
63 # Select all atoms
64 file.write('set all [atomselect top "all"]\n')
65 # Write the current topology in 'pdb' format
66 file.write('$all frame first\n')
67 file.write(f'$all writepdb {output_structure_filename}\n')
68 file.write('exit\n')
70 # Run VMD
71 vmd_process = run([
72 "vmd",
73 input_structure_filename,
74 "-e",
75 commands_filename,
76 "-dispdev",
77 "none"
78 ], stdout=PIPE, stderr=PIPE)
79 output_logs = vmd_process.stdout.decode()
81 if not exists(output_structure_filename):
82 print(output_logs)
83 error_logs = vmd_process.stderr.decode()
84 print(error_logs)
85 raise SystemExit('Something went wrong with VMD')
87 os.remove(commands_filename)
88# Set function supported formats
89vmd_to_pdb.format_sets = [
90 {
91 'inputs': {
92 'input_structure_filename': {'psf', 'parm', 'prmtop'},
93 'input_trajectory_filename': vmd_supported_trajectory_formats
94 },
95 'outputs': {
96 'output_structure_filename': {'pdb'}
97 }
98 }
99]
101# This tool allows you to set the chain of all atoms in a selection
102# This is powered by VMD and thus the selection lenguage must be the VMD's
103# Arguments are as follows:
104# 1 - Input pdb filename
105# 2 - Atom selection (All atoms by defualt)
106# 3 - Chain letter (May be the flag 'fragment', which is the default indeed)
107# 4 - Output pdb filename (Input filename by default)
108# WARNING: When no selection is passed, if only a part of a "fragment" is missing the chain then the whole fragment will be affected
109# WARNING: VMD only handles fragments if there are less fragments than letters in the alphabet
110# DEPRECATED: Use the structures chainer instead
111def chainer (
112 input_pdb_filename : str,
113 atom_selection : Optional[str] = None,
114 chain_letter : Optional[str] = None,
115 output_pdb_filename : Optional[str] = None
116):
118 # If no atom selection is provided then all atoms are selected
119 if not atom_selection:
120 atom_selection = 'all'
122 # Escape TCL meaningful characters
123 escaped_atom_selection = escape_tcl_selection(atom_selection)
125 # If no chain letter is provided then the flag 'fragment' is used
126 if not chain_letter:
127 chain_letter = 'fragment'
129 # If no output filename is provided then use input filename as output filename
130 if not output_pdb_filename:
131 output_pdb_filename = input_pdb_filename
133 # Check the file exists
134 if not exists(input_pdb_filename):
135 raise SystemExit('ERROR: The file does not exist')
137 with open(commands_filename, "w") as file:
138 # Select the specified atoms and set the specified chain
139 file.write(f'set atoms [atomselect top "{escaped_atom_selection}"]\n')
140 # In case chain letter is not a letter but the 'fragment' flag, asign chains by fragment
141 # Fragments are atom groups which are not connected by any bond
142 if chain_letter == 'fragment':
143 # Get all different chain names
144 file.write('set chains_sample [lsort -unique [${atoms} get chain]]\n')
145 # Set letters in alphabetic order
146 file.write('set letters "A B C D E F G H I J K L M N O P Q R S T U V W X Y Z"\n')
147 # Get the number of fragments
148 file.write('set fragment_number [llength [lsort -unique -integer [${atoms} get fragment]]]\n')
149 # For each fragment, set the chain of all atoms which belong to this fragment alphabetically
150 # e.g. fragment 0 -> chain A, fragment 1 -> chain B, ...
151 file.write('for {set i 0} {$i <= $fragment_number} {incr i} {\n')
152 file.write(' set fragment_atoms [atomselect top "fragment $i"]\n')
153 file.write(' $fragment_atoms set chain [lindex $letters $i]\n')
154 file.write('}\n')
155 # Otherwise, set the specified chain
156 else:
157 file.write(f'$atoms set chain {chain_letter}\n')
158 # Write the current topology in 'pdb' format
159 file.write('set all [atomselect top "all"]\n')
160 file.write('$all frame first\n')
161 file.write(f'$all writepdb {output_pdb_filename}\n')
162 file.write('exit\n')
164 # Run VMD
165 logs = run([
166 "vmd",
167 input_pdb_filename,
168 "-e",
169 commands_filename,
170 "-dispdev",
171 "none"
172 ], stdout=PIPE, stderr=PIPE).stdout.decode()
174 # If the expected output file was not generated then stop here and warn the user
175 if not exists(output_pdb_filename):
176 print(logs)
177 raise SystemExit('Something went wrong with VMD')
179 # Remove the vmd commands file
180 os.remove(commands_filename)
182# Get vmd supported trajectories merged and converted to a different format
183# WARNING: Note that this process is not memory efficient so beware the size of trajectories to be converted
184# WARNING: The input structure filename may be None
185def merge_and_convert_trajectories (
186 input_structure_filename : Optional[str],
187 input_trajectory_filenames : List[str],
188 output_trajectory_filename : str
189 ):
191 warn('You are using a not memory efficient tool. If the trajectory is too big your system may not hold it.')
192 print('Note that we cannot display the progress of the conversion since we are using VMD')
194 # Get the format to export coordinates
195 output_trajectory_file = File(output_trajectory_filename)
196 output_trajectory_format = output_trajectory_file.format
198 # Although 'crd' and 'mdcrd' may be the same, VMD only recognizes 'crd' as exporting coordinates file type
199 if output_trajectory_format == 'mdcrd':
200 output_trajectory_format = 'crd'
202 # Prepare a script for the VMD to automate the data parsing. This is Tcl lenguage
203 # In addition, if chains are missing, this script asigns chains by fragment
204 # Fragments are atom groups which are not connected by any bond
205 with open(commands_filename, "w") as file:
206 # Select all atoms
207 file.write('set all [atomselect top "all"]\n')
208 # Write the current trajectory in the specified format format
209 file.write(f'animate write {output_trajectory_format} {output_trajectory_filename} waitfor all sel $all\n')
210 file.write('exit\n')
212 inputs = [ input_structure_filename, *input_trajectory_filenames ] if input_structure_filename else input_trajectory_filenames
214 # Run VMD
215 logs = run([
216 "vmd",
217 *inputs,
218 "-e",
219 commands_filename,
220 "-dispdev",
221 "none"
222 ], stdout=PIPE, stderr=STDOUT).stdout.decode() # Redirect errors to the output in order to dont show them in console
224 # If the expected output file was not generated then stop here and warn the user
225 if not exists(output_trajectory_filename):
226 print(logs)
227 raise SystemExit('Something went wrong with VMD')
229 # Remove the vmd commands file
230 os.remove(commands_filename)
232# Set function supported formats
233merge_and_convert_trajectories.format_sets = [
234 {
235 'inputs': {
236 'input_structure_filename': None,
237 'input_trajectory_filenames': {'dcd', 'xtc', 'trr', 'nc'}
238 },
239 'outputs': {
240 # NEVER FORGET: VMD cannot write to all formats it supports to read
241 'output_trajectory_filename': {'mdcrd', 'crd', 'dcd', 'trr'}
242 }
243 },
244 {
245 'inputs': {
246 'input_structure_filename': vmd_supported_structure_formats,
247 'input_trajectory_filenames': {'mdcrd', 'crd', 'rst7'}
248 },
249 'outputs': {
250 # NEVER FORGET: VMD cannot write to all formats it supports to read
251 'output_trajectory_filename': {'mdcrd', 'crd', 'dcd', 'trr'}
252 }
253 }
254]
256# Given an atom selection in vmd syntax, return the list of atom indices it corresponds to
257def get_vmd_selection_atom_indices (input_structure_filename : str, selection : str) -> List[int]:
259 # Escape TCL meaningful characters
260 escaped_selection = escape_tcl_selection(selection)
262 # Prepare a script for VMD to run. This is Tcl language
263 # The output of the script will be written to a txt file
264 atom_indices_filename = '.vmd_output.txt'
265 with open(commands_filename, "w") as file:
266 # Select the specified atoms
267 file.write(f'set selection [atomselect top "{escaped_selection}"]\n')
268 # Save atom indices from the selection
269 file.write('set indices [$selection list]\n')
270 # Write atom indices to a file
271 file.write(f'set indices_file [open {atom_indices_filename} w]\n')
272 file.write('puts $indices_file $indices\n')
273 file.write('exit\n')
275 # Run VMD
276 logs = run([
277 "vmd",
278 input_structure_filename,
279 "-e",
280 commands_filename,
281 "-dispdev",
282 "none"
283 ], stdout=PIPE, stderr=PIPE).stdout.decode()
285 # If the expected output file was not generated then stop here and warn the user
286 if not exists(atom_indices_filename):
287 print(logs)
288 raise SystemExit('Something went wrong with VMD')
290 # Read the VMD output
291 with open(atom_indices_filename, 'r') as file:
292 raw_atom_indices = file.read()
294 # Parse the atom indices string to an array of integers
295 atom_indices = [ int(i) for i in raw_atom_indices.split() ]
297 # Remove trahs files
298 trash_files = [ commands_filename, atom_indices_filename ]
299 for trash_file in trash_files:
300 os.remove(trash_file)
302 return atom_indices
304# Set a function to retrieve all covalent (strong) bonds in a structure
305# You may provide an atom selection as well
306def get_covalent_bonds (structure_filename : str, selection : Optional['Selection'] = None) -> List[ List[int] ]:
308 # Parse the selection to vmd
309 vmd_selection = 'all'
310 if selection:
311 vmd_selection = selection.to_vmd()
313 # Prepare a script for the VMD to automate the commands. This is Tcl lenguage
314 output_bonds_file = '.bonds.txt'
315 with open(commands_filename, "w") as file:
316 # Select atoms
317 file.write(f'set atoms [atomselect top "{vmd_selection}"]\n')
318 # Save covalent bonds
319 file.write('set bonds [$atoms getbonds]\n')
320 # Write those bonds to a file
321 file.write(f'set bondsfile [open {output_bonds_file} w]\n')
322 file.write('puts $bondsfile $bonds\n')
323 file.write('exit\n')
325 # Run VMD
326 logs = run([
327 "vmd",
328 structure_filename,
329 "-e",
330 commands_filename,
331 "-dispdev",
332 "none"
333 ], stdout=PIPE, stderr=PIPE).stdout.decode()
335 # If the output file is missing at this point then it means something went wrong
336 if not exists(output_bonds_file):
337 print(logs)
338 raise SystemExit('Something went wrong with VMD')
340 # Read the VMD output
341 with open(output_bonds_file, 'r') as file:
342 raw_bonds = file.read()
344 # Remove vmd files since they are no longer usefull
345 for f in [ commands_filename, output_bonds_file ]:
346 os.remove(f)
348 # Sometimes there is a breakline at the end of the raw bonds string and it must be removed
349 # Add a space at the end of the string to make the parser get the last character
350 raw_bonds = raw_bonds.replace('\n', '') + ' '
352 # Parse the raw bonds string to a list of atom bonds (i.e. a list of lists of integers)
353 # Raw bonds format is (for each atom in the selection):
354 # '{index1, index2, index3 ...}' with the index of each connected atom
355 # 'index' if there is only one connected atom
356 # '{}' if there are no connected atoms
357 bonds_per_atom = []
358 last_atom_index = ''
359 last_atom_bonds = []
360 in_brackets = False
361 for character in raw_bonds:
362 if character == ' ':
363 if len(last_atom_index) > 0:
364 if in_brackets:
365 last_atom_bonds.append(int(last_atom_index))
366 else:
367 bonds_per_atom.append([int(last_atom_index)])
368 last_atom_index = ''
369 continue
370 if character == '{':
371 in_brackets = True
372 continue
373 if character == '}':
374 if last_atom_index == '':
375 bonds_per_atom.append([])
376 in_brackets = False
377 continue
378 last_atom_bonds.append(int(last_atom_index))
379 last_atom_index = ''
380 bonds_per_atom.append(last_atom_bonds)
381 last_atom_bonds = []
382 in_brackets = False
383 continue
384 last_atom_index += character
386 return bonds_per_atom
388# Set a function to retrieve covalent (strong) bonds between 2 atom selections
389def get_covalent_bonds_between (
390 structure_filename : str,
391 # Selections may be either selection instances or selection strings already in VMD format
392 selection_1 : Union['Selection', str],
393 selection_2 : Union['Selection', str]
394) -> List[ List[int] ]:
396 # Parse selections (if not parsed yet)
397 parsed_selection_1 = selection_1 if type(selection_1) == str else selection_1.to_vmd()
398 parsed_selection_2 = selection_2 if type(selection_2) == str else selection_2.to_vmd()
400 # Prepare a script for the VMD to automate the commands. This is Tcl lenguage
401 output_index_1_file = '.index1.txt'
402 output_index_2_file = '.index2.txt'
403 output_bonds_file = '.bonds.ext'
404 with open(commands_filename, "w") as file:
405 # Select the specified atoms in selection 1
406 file.write(f'set sel1 [atomselect top "{parsed_selection_1}"]\n')
407 # Save all atom index in the selection
408 file.write('set index1 [$sel1 list]\n')
409 # Write those index to a file
410 file.write(f'set indexfile1 [open {output_index_1_file} w]\n')
411 file.write('puts $indexfile1 $index1\n')
412 # Save all covalent bonds in the selection
413 file.write('set bonds [$sel1 getbonds]\n')
414 # Write those bonds to a file
415 file.write(f'set bondsfile [open {output_bonds_file} w]\n')
416 file.write('puts $bondsfile $bonds\n')
417 # Select the specified atoms in selection 2
418 file.write(f'set sel2 [atomselect top "{parsed_selection_2}"]\n')
419 # Save all atom index in the selection
420 file.write('set index2 [$sel2 list]\n')
421 # Write those index to a file
422 file.write(f'set indexfile2 [open {output_index_2_file} w]\n')
423 file.write('puts $indexfile2 $index2\n')
424 file.write('exit\n')
426 # Run VMD
427 logs = run([
428 "vmd",
429 structure_filename,
430 "-e",
431 commands_filename,
432 "-dispdev",
433 "none"
434 ], stdout=PIPE, stderr=PIPE).stdout.decode()
436 # If the expected output file was not generated then stop here and warn the user
437 if not exists(output_index_1_file) or not exists(output_bonds_file) or not exists(output_index_2_file):
438 print(logs)
439 raise SystemExit('Something went wrong with VMD')
441 # Read the VMD output
442 with open(output_index_1_file, 'r') as file:
443 raw_index_1 = file.read()
444 with open(output_bonds_file, 'r') as file:
445 raw_bonds = file.read()
446 with open(output_index_2_file, 'r') as file:
447 raw_index_2 = file.read()
449 # Remove vmd files since they are no longer usefull
450 for f in [ commands_filename, output_index_1_file, output_index_2_file, output_bonds_file ]:
451 os.remove(f)
453 # Sometimes there is a breakline at the end of the raw bonds string and it must be removed
454 # Add a space at the end of the string to make the parser get the last character
455 raw_bonds = raw_bonds.replace('\n', '') + ' '
457 # Raw indexes is a string with all indexes separated by spaces
458 index_1 = [ int(i) for i in raw_index_1.split() ]
459 index_2 = [ int(i) for i in raw_index_2.split() ]
461 # Parse the raw bonds string to a list of atom bonds (i.e. a list of lists of integers)
462 # Raw bonds format is (for each atom in the selection):
463 # '{index1, index2, index3 ...}' with the index of each connected atom
464 # 'index' if there is only one connected atom
465 # '{}' if there are no connected atoms
466 bonds_per_atom = []
467 last_atom_index = ''
468 last_atom_bonds = []
469 in_brackets = False
470 for character in raw_bonds:
471 if character == ' ':
472 if len(last_atom_index) > 0:
473 if in_brackets:
474 last_atom_bonds.append(int(last_atom_index))
475 else:
476 bonds_per_atom.append([int(last_atom_index)])
477 last_atom_index = ''
478 continue
479 if character == '{':
480 in_brackets = True
481 continue
482 if character == '}':
483 if last_atom_index == '':
484 bonds_per_atom.append([])
485 in_brackets = False
486 continue
487 last_atom_bonds.append(int(last_atom_index))
488 last_atom_index = ''
489 bonds_per_atom.append(last_atom_bonds)
490 last_atom_bonds = []
491 in_brackets = False
492 continue
493 last_atom_index += character
495 # At this point indexes and bonds from the first selection should match in number
496 if len(index_1) != len(bonds_per_atom):
497 raise ValueError(f'Indexes ({len(index_1)}) and atom bonds ({len(bonds_per_atom)}) do not match in number')
499 # Now get all covalent bonds which include an index from the atom selection 2
500 crossed_bonds = []
501 for i, index in enumerate(index_1):
502 bonds = bonds_per_atom[i]
503 for bond in bonds:
504 if bond in index_2:
505 crossed_bond = (index, bond)
506 crossed_bonds.append(crossed_bond)
508 return crossed_bonds
510# Given two atom selections, find interface atoms and return their indices
511# Interface atoms are those atoms closer than the cutoff in at least 1 frame along a trajectory
512# Return also atom indices for the whole selections
513def get_interface_atom_indices (
514 input_structure_filepath : str,
515 input_trajectory_filepath : str,
516 selection_1 : str,
517 selection_2 : str,
518 distance_cutoff : float,
519) -> List[int]:
521 # Set the interface selections
522 interface_selection_1 = (f'({selection_1}) and within {distance_cutoff} of ({selection_2})')
523 interface_selection_2 = (f'({selection_2 }) and within {distance_cutoff} of ({selection_1})')
525 # Set the output txt files for vmd to write the atom indices
526 # Note that these output files are deleted at the end of this function
527 selection_1_filename = '.selection_1.txt'
528 selection_2_filename = '.selection_2.txt'
529 interface_selection_1_filename = '.interface_selection_1.txt'
530 interface_selection_2_filename = '.interface_selection_2.txt'
531 interacting_frames_filename = '.iframes.txt'
532 total_frames_filename = '.nframes.txt'
534 # Prepare a script for VMD to run. This is Tcl language
535 commands_filename = '.commands.vmd'
536 with open(commands_filename, "w") as file:
537 # -------------------------------------------
538 # First get the whole selection atom indices
539 # -------------------------------------------
540 # Select the specified atoms
541 file.write(f'set selection [atomselect top "{selection_1}"]\n')
542 # Save atom indices from the selection
543 file.write('set indices [$selection list]\n')
544 # Write atom indices to a file
545 file.write(f'set indices_file [open {selection_1_filename} w]\n')
546 file.write('puts $indices_file $indices\n')
547 # Select the specified atoms
548 file.write(f'set selection [atomselect top "{selection_2}"]\n')
549 # Save atom indices from the selection
550 file.write('set indices [$selection list]\n')
551 # Write atom indices to a file
552 file.write(f'set indices_file [open {selection_2_filename} w]\n')
553 file.write('puts $indices_file $indices\n')
554 # -------------------------------------------
555 # Now get the interface selection atom indices
556 # Also count the number of frames where there is at least one interacting residue
557 # -------------------------------------------
558 # Capture indices for each frame in the trajectory
559 file.write('set accumulated_interface1_atom_indices []\n')
560 file.write('set accumulated_interface2_atom_indices []\n')
561 file.write(f'set interface1 [atomselect top "{interface_selection_1}"]\n')
562 file.write(f'set interface2 [atomselect top "{interface_selection_2}"]\n')
563 # Capture the number of frames where the interaction happens
564 file.write('set iframes 0\n')
565 # Get the number of frames in the trajectory
566 file.write('set nframes [molinfo top get numframes]\n')
567 # Iterate over each frame
568 # Note that we skip the first frame (i = 1, not i = 0) since it belongs to the structure
569 file.write('for { set i 1 } { $i < $nframes } { incr i } {\n')
570 # Update the selection in the current frame
571 file.write(' $interface1 frame $i\n')
572 file.write(' $interface1 update\n')
573 # Add its atom indices to the acumulated atom indices
574 file.write(' set interface1_atom_indices [$interface1 list]\n')
575 file.write(' set accumulated_interface1_atom_indices [concat $accumulated_interface1_atom_indices $interface1_atom_indices ]\n')
576 # Repeat with the selection 2
577 file.write(' $interface2 frame $i\n')
578 file.write(' $interface2 update\n')
579 file.write(' set interface2_atom_indices [$interface2 list]\n')
580 file.write(' set accumulated_interface2_atom_indices [concat $accumulated_interface2_atom_indices $interface2_atom_indices ]\n')
581 # If there was at least one atom in one of the interactions then add one to the interaction frame count
582 # Note that checking both interactions would be redundant so one is enough
583 file.write(' if { [llength $interface1_atom_indices] > 0 } {\n')
584 file.write(' incr iframes\n')
585 file.write(' }\n')
586 file.write('}\n')
587 # Write the number of interacting frames and total frames to files
588 file.write(f'set iframes_file [open {interacting_frames_filename} w]\n')
589 file.write('puts $iframes_file $iframes\n')
590 file.write(f'set nframes_file [open {total_frames_filename} w]\n')
591 # Note that we must substract 1 from the total frames count since the first frame was skipped
592 file.write('set cframes [expr $nframes - 1]\n')
593 file.write('puts $nframes_file $cframes\n')
594 # Remove duplicated indices
595 # Use the -integer argument to do the correct sorting
596 # Otherwise numbers are sorted as strings so you have 1, 10, 100, 2, etc.
597 file.write('set unique_interface1_atom_indices [lsort -integer -unique $accumulated_interface1_atom_indices]\n')
598 file.write('set unique_interface2_atom_indices [lsort -integer -unique $accumulated_interface2_atom_indices]\n')
599 # Write indices to files
600 file.write(f'set indices_file [open {interface_selection_1_filename} w]\n')
601 file.write('puts $indices_file $unique_interface1_atom_indices\n')
602 file.write(f'set indices_file [open {interface_selection_2_filename} w]\n')
603 file.write('puts $indices_file $unique_interface2_atom_indices\n')
604 file.write('exit\n')
606 # Run VMD
607 logs = run([
608 "vmd",
609 input_structure_filepath,
610 input_trajectory_filepath,
611 "-e",
612 commands_filename,
613 "-dispdev",
614 "none"
615 ], stdout=PIPE, stderr=PIPE).stdout.decode()
617 # If any of the output files do not exist at this point then it means something went wrong with vmd
618 expected_output_files = [
619 selection_1_filename,
620 selection_2_filename,
621 interface_selection_1_filename,
622 interface_selection_2_filename,
623 interacting_frames_filename,
624 total_frames_filename
625 ]
626 for output_file in expected_output_files:
627 if not os.path.exists(output_file):
628 print(logs)
629 raise SystemExit('Something went wrong with VMD')
631 # Set a function to read the VMD output and parse the atom indices string to an array of integers
632 def process_vmd_output (output_filename : str) -> List[int]:
633 with open(output_filename, 'r') as file:
634 raw_atom_indices = file.read()
635 return [ int(i) for i in raw_atom_indices.split() ]
637 # Read the VMD output
638 selection_1_atom_indices = process_vmd_output(selection_1_filename)
639 selection_2_atom_indices = process_vmd_output(selection_2_filename)
640 selection_1_interface_atom_indices = process_vmd_output(interface_selection_1_filename)
641 selection_2_interface_atom_indices = process_vmd_output(interface_selection_2_filename)
642 interacting_frames = process_vmd_output(interacting_frames_filename)[0]
643 total_frames = process_vmd_output(total_frames_filename)[0]
645 # Remove trash files
646 trash_files = [ commands_filename ] + expected_output_files
647 for trash_file in trash_files:
648 os.remove(trash_file)
650 # Return the results
651 return {
652 'selection_1_atom_indices': selection_1_atom_indices,
653 'selection_2_atom_indices': selection_2_atom_indices,
654 'selection_1_interface_atom_indices': selection_1_interface_atom_indices,
655 'selection_2_interface_atom_indices': selection_2_interface_atom_indices,
656 'interacting_frames': interacting_frames,
657 'total_frames': total_frames
658 }