Coverage for mddb_workflow/utils/vmd_spells.py: 79%
292 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-10-29 15:48 +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 mddb_workflow.utils.file import File
11from mddb_workflow.utils.type_hints import *
12from mddb_workflow.utils.auxiliar import warn
14# Set characters to be escaped since they have a meaning in TCL
15TCL_RESERVED_CHARACTERS = ['"','[',']']
17def escape_tcl_selection (selection : str) -> str:
18 """ Given a VMD atom selection string, escape TCL meaningful characters and return the escaped string. """
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]
101def chainer (
102 input_pdb_filename : str,
103 atom_selection : Optional[str] = None,
104 chain_letter : Optional[str] = None,
105 output_pdb_filename : Optional[str] = None
106):
107 """ This tool allows you to set the chain of all atoms in a selection.
108 This is powered by VMD and thus the selection lenguage must be the VMD's.
109 Arguments are as follows:
111 1. Input pdb filename
112 2. Atom selection (All atoms by defualt)
113 3. Chain letter (May be the flag 'fragment', which is the default indeed)
114 4. Output pdb filename (Input filename by default)
116 WARNING: When no selection is passed, if only a part of a "fragment" is missing the chain then the whole fragment will be affected
117 WARNING: VMD only handles fragments if there are less fragments than letters in the alphabet
118 DEPRECATED: Use the structures chainer instead. """
120 # If no atom selection is provided then all atoms are selected
121 if not atom_selection:
122 atom_selection = 'all'
124 # Escape TCL meaningful characters
125 escaped_atom_selection = escape_tcl_selection(atom_selection)
127 # If no chain letter is provided then the flag 'fragment' is used
128 if not chain_letter:
129 chain_letter = 'fragment'
131 # If no output filename is provided then use input filename as output filename
132 if not output_pdb_filename:
133 output_pdb_filename = input_pdb_filename
135 # Check the file exists
136 if not exists(input_pdb_filename):
137 raise SystemExit('ERROR: The file does not exist')
139 with open(commands_filename, "w") as file:
140 # Select the specified atoms and set the specified chain
141 file.write(f'set atoms [atomselect top "{escaped_atom_selection}"]\n')
142 # In case chain letter is not a letter but the 'fragment' flag, asign chains by fragment
143 # Fragments are atom groups which are not connected by any bond
144 if chain_letter == 'fragment':
145 # Get all different chain names
146 file.write('set chains_sample [lsort -unique [${atoms} get chain]]\n')
147 # Set letters in alphabetic order
148 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')
149 # Get the number of fragments
150 file.write('set fragment_number [llength [lsort -unique -integer [${atoms} get fragment]]]\n')
151 # For each fragment, set the chain of all atoms which belong to this fragment alphabetically
152 # e.g. fragment 0 -> chain A, fragment 1 -> chain B, ...
153 file.write('for {set i 0} {$i <= $fragment_number} {incr i} {\n')
154 file.write(' set fragment_atoms [atomselect top "fragment $i"]\n')
155 file.write(' $fragment_atoms set chain [lindex $letters $i]\n')
156 file.write('}\n')
157 # Otherwise, set the specified chain
158 else:
159 file.write(f'$atoms set chain {chain_letter}\n')
160 # Write the current topology in 'pdb' format
161 file.write('set all [atomselect top "all"]\n')
162 file.write('$all frame first\n')
163 file.write(f'$all writepdb {output_pdb_filename}\n')
164 file.write('exit\n')
166 # Run VMD
167 logs = run([
168 "vmd",
169 input_pdb_filename,
170 "-e",
171 commands_filename,
172 "-dispdev",
173 "none"
174 ], stdout=PIPE, stderr=PIPE).stdout.decode()
176 # If the expected output file was not generated then stop here and warn the user
177 if not exists(output_pdb_filename):
178 print(logs)
179 raise SystemExit('Something went wrong with VMD')
181 # Remove the vmd commands file
182 os.remove(commands_filename)
184def merge_and_convert_trajectories (
185 input_structure_filename : Optional[str],
186 input_trajectory_filenames : list[str],
187 output_trajectory_filename : str
188 ):
189 """ Get vmd supported trajectories merged and converted to a different format.
190 WARNING: Note that this process is not memory efficient so beware the size of trajectories to be converted.
191 WARNING: The input structure filename may be None. """
193 warn('You are using a not memory efficient tool. If the trajectory is too big your system may not hold it.')
194 print('Note that we cannot display the progress of the conversion since we are using VMD')
196 # Get the format to export coordinates
197 output_trajectory_file = File(output_trajectory_filename)
198 output_trajectory_format = output_trajectory_file.format
200 # Although 'crd' and 'mdcrd' may be the same, VMD only recognizes 'crd' as exporting coordinates file type
201 if output_trajectory_format == 'mdcrd':
202 output_trajectory_format = 'crd'
204 # Prepare a script for the VMD to automate the data parsing. This is Tcl lenguage
205 # In addition, if chains are missing, this script asigns chains by fragment
206 # Fragments are atom groups which are not connected by any bond
207 with open(commands_filename, "w") as file:
208 # Select all atoms
209 file.write('set all [atomselect top "all"]\n')
210 # Write the current trajectory in the specified format format
211 file.write(f'animate write {output_trajectory_format} {output_trajectory_filename} waitfor all sel $all\n')
212 file.write('exit\n')
214 inputs = [ input_structure_filename, *input_trajectory_filenames ] if input_structure_filename else input_trajectory_filenames
216 # Run VMD
217 logs = run([
218 "vmd",
219 *inputs,
220 "-e",
221 commands_filename,
222 "-dispdev",
223 "none"
224 ], stdout=PIPE, stderr=STDOUT).stdout.decode() # Redirect errors to the output in order to dont show them in console
226 # If the expected output file was not generated then stop here and warn the user
227 if not exists(output_trajectory_filename):
228 print(logs)
229 raise SystemExit('Something went wrong with VMD')
231 # Remove the vmd commands file
232 os.remove(commands_filename)
234# Set function supported formats
235merge_and_convert_trajectories.format_sets = [
236 {
237 'inputs': {
238 'input_structure_filename': None,
239 'input_trajectory_filenames': {'dcd', 'xtc', 'trr', 'nc'}
240 },
241 'outputs': {
242 # NEVER FORGET: VMD cannot write to all formats it supports to read
243 'output_trajectory_filename': {'mdcrd', 'crd', 'dcd', 'trr'}
244 }
245 },
246 {
247 'inputs': {
248 'input_structure_filename': vmd_supported_structure_formats,
249 'input_trajectory_filenames': {'mdcrd', 'crd', 'rst7'}
250 },
251 'outputs': {
252 # NEVER FORGET: VMD cannot write to all formats it supports to read
253 'output_trajectory_filename': {'mdcrd', 'crd', 'dcd', 'trr'}
254 }
255 }
256]
258def get_vmd_selection_atom_indices (input_structure_filename : str, selection : str) -> list[int]:
259 """ Given an atom selection in vmd syntax, return the list of atom indices it corresponds to. """
261 # Escape TCL meaningful characters
262 escaped_selection = escape_tcl_selection(selection)
264 # Prepare a script for VMD to run. This is Tcl language
265 # The output of the script will be written to a txt file
266 atom_indices_filename = '.vmd_output.txt'
267 with open(commands_filename, "w") as file:
268 # Select the specified atoms
269 file.write(f'set selection [atomselect top "{escaped_selection}"]\n')
270 # Save atom indices from the selection
271 file.write('set indices [$selection list]\n')
272 # Write atom indices to a file
273 file.write(f'set indices_file [open {atom_indices_filename} w]\n')
274 file.write('puts $indices_file $indices\n')
275 file.write('exit\n')
277 # Run VMD
278 logs = run([
279 "vmd",
280 input_structure_filename,
281 "-e",
282 commands_filename,
283 "-dispdev",
284 "none"
285 ], stdout=PIPE, stderr=PIPE).stdout.decode()
287 # If the expected output file was not generated then stop here and warn the user
288 if not exists(atom_indices_filename):
289 print(logs)
290 raise SystemExit('Something went wrong with VMD')
292 # Read the VMD output
293 with open(atom_indices_filename, 'r') as file:
294 raw_atom_indices = file.read()
296 # Parse the atom indices string to an array of integers
297 atom_indices = [ int(i) for i in raw_atom_indices.split() ]
299 # Remove trahs files
300 trash_files = [ commands_filename, atom_indices_filename ]
301 for trash_file in trash_files:
302 os.remove(trash_file)
304 return atom_indices
306def get_covalent_bonds (structure_filename : str, selection : Optional['Selection'] = None) -> list[ list[int] ]:
307 """ Set a function to retrieve all covalent (strong) bonds in a structure using VMD.
308 You may provide an atom selection as well. """
310 # Parse the selection to vmd
311 vmd_selection = 'all'
312 if selection:
313 vmd_selection = selection.to_vmd()
315 # Prepare a script for the VMD to automate the commands. This is Tcl lenguage
316 output_bonds_file = '.bonds.txt'
317 with open(commands_filename, "w") as file:
318 # Select atoms
319 file.write(f'set atoms [atomselect top "{vmd_selection}"]\n')
320 # Save covalent bonds
321 file.write('set bonds [$atoms getbonds]\n')
322 # Write those bonds to a file
323 file.write(f'set bondsfile [open {output_bonds_file} w]\n')
324 file.write('puts $bondsfile $bonds\n')
325 file.write('exit\n')
327 # Run VMD
328 logs = run([
329 "vmd",
330 structure_filename,
331 "-e",
332 commands_filename,
333 "-dispdev",
334 "none"
335 ], stdout=PIPE, stderr=PIPE).stdout.decode()
337 # If the output file is missing at this point then it means something went wrong
338 if not exists(output_bonds_file):
339 print(logs)
340 raise SystemExit('Something went wrong with VMD')
342 # Read the VMD output
343 with open(output_bonds_file, 'r') as file:
344 raw_bonds = file.read()
346 # Remove vmd files since they are no longer usefull
347 for f in [ commands_filename, output_bonds_file ]:
348 os.remove(f)
350 # Sometimes there is a breakline at the end of the raw bonds string and it must be removed
351 # Add a space at the end of the string to make the parser get the last character
352 raw_bonds = raw_bonds.replace('\n', '') + ' '
354 # Parse the raw bonds string to a list of atom bonds (i.e. a list of lists of integers)
355 # Raw bonds format is (for each atom in the selection):
356 # '{index1, index2, index3 ...}' with the index of each connected atom
357 # 'index' if there is only one connected atom
358 # '{}' if there are no connected atoms
359 bonds_per_atom = []
360 last_atom_index = ''
361 last_atom_bonds = []
362 in_brackets = False
363 for character in raw_bonds:
364 if character == ' ':
365 if len(last_atom_index) > 0:
366 if in_brackets:
367 last_atom_bonds.append(int(last_atom_index))
368 else:
369 bonds_per_atom.append([int(last_atom_index)])
370 last_atom_index = ''
371 continue
372 if character == '{':
373 in_brackets = True
374 continue
375 if character == '}':
376 if last_atom_index == '':
377 bonds_per_atom.append([])
378 in_brackets = False
379 continue
380 last_atom_bonds.append(int(last_atom_index))
381 last_atom_index = ''
382 bonds_per_atom.append(last_atom_bonds)
383 last_atom_bonds = []
384 in_brackets = False
385 continue
386 last_atom_index += character
388 return bonds_per_atom
390def get_covalent_bonds_between (
391 structure_filename : str,
392 # Selections may be either selection instances or selection strings already in VMD format
393 selection_1 : Union['Selection', str],
394 selection_2 : Union['Selection', str]
395) -> list[ list[int] ]:
396 """ Set a function to retrieve covalent (strong) bonds between 2 atom selections. """
398 # Parse selections (if not parsed yet)
399 parsed_selection_1 = selection_1 if type(selection_1) == str else selection_1.to_vmd()
400 parsed_selection_2 = selection_2 if type(selection_2) == str else selection_2.to_vmd()
402 # Prepare a script for the VMD to automate the commands. This is Tcl lenguage
403 output_index_1_file = '.index1.txt'
404 output_index_2_file = '.index2.txt'
405 output_bonds_file = '.bonds.ext'
406 with open(commands_filename, "w") as file:
407 # Select the specified atoms in selection 1
408 file.write(f'set sel1 [atomselect top "{parsed_selection_1}"]\n')
409 # Save all atom index in the selection
410 file.write('set index1 [$sel1 list]\n')
411 # Write those index to a file
412 file.write(f'set indexfile1 [open {output_index_1_file} w]\n')
413 file.write('puts $indexfile1 $index1\n')
414 # Save all covalent bonds in the selection
415 file.write('set bonds [$sel1 getbonds]\n')
416 # Write those bonds to a file
417 file.write(f'set bondsfile [open {output_bonds_file} w]\n')
418 file.write('puts $bondsfile $bonds\n')
419 # Select the specified atoms in selection 2
420 file.write(f'set sel2 [atomselect top "{parsed_selection_2}"]\n')
421 # Save all atom index in the selection
422 file.write('set index2 [$sel2 list]\n')
423 # Write those index to a file
424 file.write(f'set indexfile2 [open {output_index_2_file} w]\n')
425 file.write('puts $indexfile2 $index2\n')
426 file.write('exit\n')
428 # Run VMD
429 logs = run([
430 "vmd",
431 structure_filename,
432 "-e",
433 commands_filename,
434 "-dispdev",
435 "none"
436 ], stdout=PIPE, stderr=PIPE).stdout.decode()
438 # If the expected output file was not generated then stop here and warn the user
439 if not exists(output_index_1_file) or not exists(output_bonds_file) or not exists(output_index_2_file):
440 print(logs)
441 raise SystemExit('Something went wrong with VMD')
443 # Read the VMD output
444 with open(output_index_1_file, 'r') as file:
445 raw_index_1 = file.read()
446 with open(output_bonds_file, 'r') as file:
447 raw_bonds = file.read()
448 with open(output_index_2_file, 'r') as file:
449 raw_index_2 = file.read()
451 # Remove vmd files since they are no longer usefull
452 for f in [ commands_filename, output_index_1_file, output_index_2_file, output_bonds_file ]:
453 os.remove(f)
455 # Sometimes there is a breakline at the end of the raw bonds string and it must be removed
456 # Add a space at the end of the string to make the parser get the last character
457 raw_bonds = raw_bonds.replace('\n', '') + ' '
459 # Raw indexes is a string with all indexes separated by spaces
460 index_1 = [ int(i) for i in raw_index_1.split() ]
461 index_2 = set([ int(i) for i in raw_index_2.split() ])
463 # Parse the raw bonds string to a list of atom bonds (i.e. a list of lists of integers)
464 # Raw bonds format is (for each atom in the selection):
465 # '{index1, index2, index3 ...}' with the index of each connected atom
466 # 'index' if there is only one connected atom
467 # '{}' if there are no connected atoms
468 bonds_per_atom = []
469 last_atom_index = ''
470 last_atom_bonds = []
471 in_brackets = False
472 for character in raw_bonds:
473 if character == ' ':
474 if len(last_atom_index) > 0:
475 if in_brackets:
476 last_atom_bonds.append(int(last_atom_index))
477 else:
478 bonds_per_atom.append([int(last_atom_index)])
479 last_atom_index = ''
480 continue
481 if character == '{':
482 in_brackets = True
483 continue
484 if character == '}':
485 if last_atom_index == '':
486 bonds_per_atom.append([])
487 in_brackets = False
488 continue
489 last_atom_bonds.append(int(last_atom_index))
490 last_atom_index = ''
491 bonds_per_atom.append(last_atom_bonds)
492 last_atom_bonds = []
493 in_brackets = False
494 continue
495 last_atom_index += character
497 # At this point indexes and bonds from the first selection should match in number
498 if len(index_1) != len(bonds_per_atom):
499 raise ValueError(f'Indexes ({len(index_1)}) and atom bonds ({len(bonds_per_atom)}) do not match in number')
501 # Now get all covalent bonds which include an index from the atom selection 2
502 crossed_bonds = []
503 for i, index in enumerate(index_1):
504 bonds = bonds_per_atom[i]
505 for bond in bonds:
506 if bond in index_2:
507 # DANI: A tuple may make more sense than a list to define a bond
508 # DANI: However tuples are converted to lists when JSON serialized
509 # DANI: And the interactions are saved to json, so a list keeps things more coherent
510 crossed_bond = [index, bond]
511 crossed_bonds.append(crossed_bond)
513 return crossed_bonds
515def get_interface_atom_indices (
516 input_structure_filepath : str,
517 input_trajectory_filepath : str,
518 selection_1 : str,
519 selection_2 : str,
520 distance_cutoff : float,
521) -> list[int]:
522 """ Given two atom selections, find interface atoms and return their indices
523 Interface atoms are those atoms closer than the cutoff in at least 1 frame along a trajectory
524 Return also atom indices for the whole selections. """
526 # Set the interface selections
527 interface_selection_1 = (f'({selection_1}) and within {distance_cutoff} of ({selection_2})')
528 interface_selection_2 = (f'({selection_2 }) and within {distance_cutoff} of ({selection_1})')
530 # Set the output txt files for vmd to write the atom indices
531 # Note that these output files are deleted at the end of this function
532 selection_1_filename = '.selection_1.txt'
533 selection_2_filename = '.selection_2.txt'
534 interface_selection_1_filename = '.interface_selection_1.txt'
535 interface_selection_2_filename = '.interface_selection_2.txt'
536 interacting_frames_filename = '.iframes.txt'
537 total_frames_filename = '.nframes.txt'
539 # Prepare a script for VMD to run. This is Tcl language
540 commands_filename = '.commands.vmd'
541 with open(commands_filename, "w") as file:
542 # -------------------------------------------
543 # First get the whole selection atom indices
544 # -------------------------------------------
545 # Select the specified atoms
546 file.write(f'set selection [atomselect top "{selection_1}"]\n')
547 # Save atom indices from the selection
548 file.write('set indices [$selection list]\n')
549 # Write atom indices to a file
550 file.write(f'set indices_file [open {selection_1_filename} w]\n')
551 file.write('puts $indices_file $indices\n')
552 # Select the specified atoms
553 file.write(f'set selection [atomselect top "{selection_2}"]\n')
554 # Save atom indices from the selection
555 file.write('set indices [$selection list]\n')
556 # Write atom indices to a file
557 file.write(f'set indices_file [open {selection_2_filename} w]\n')
558 file.write('puts $indices_file $indices\n')
559 # -------------------------------------------
560 # Now get the interface selection atom indices
561 # Also count the number of frames where there is at least one interacting residue
562 # -------------------------------------------
563 # Capture indices for each frame in the trajectory
564 file.write('set accumulated_interface1_atom_indices []\n')
565 file.write('set accumulated_interface2_atom_indices []\n')
566 file.write(f'set interface1 [atomselect top "{interface_selection_1}"]\n')
567 file.write(f'set interface2 [atomselect top "{interface_selection_2}"]\n')
568 # Capture the number of frames where the interaction happens
569 file.write('set iframes 0\n')
570 # Get the number of frames in the trajectory
571 file.write('set nframes [molinfo top get numframes]\n')
572 # Iterate over each frame
573 # Note that we skip the first frame (i = 1, not i = 0) since it belongs to the structure
574 file.write('for { set i 1 } { $i < $nframes } { incr i } {\n')
575 # Update the selection in the current frame
576 file.write(' $interface1 frame $i\n')
577 file.write(' $interface1 update\n')
578 # Add its atom indices to the acumulated atom indices
579 file.write(' set interface1_atom_indices [$interface1 list]\n')
580 file.write(' set accumulated_interface1_atom_indices [concat $accumulated_interface1_atom_indices $interface1_atom_indices ]\n')
581 # Repeat with the selection 2
582 file.write(' $interface2 frame $i\n')
583 file.write(' $interface2 update\n')
584 file.write(' set interface2_atom_indices [$interface2 list]\n')
585 file.write(' set accumulated_interface2_atom_indices [concat $accumulated_interface2_atom_indices $interface2_atom_indices ]\n')
586 # If there was at least one atom in one of the interactions then add one to the interaction frame count
587 # Note that checking both interactions would be redundant so one is enough
588 file.write(' if { [llength $interface1_atom_indices] > 0 } {\n')
589 file.write(' incr iframes\n')
590 file.write(' }\n')
591 file.write('}\n')
592 # Write the number of interacting frames and total frames to files
593 file.write(f'set iframes_file [open {interacting_frames_filename} w]\n')
594 file.write('puts $iframes_file $iframes\n')
595 file.write(f'set nframes_file [open {total_frames_filename} w]\n')
596 # Note that we must substract 1 from the total frames count since the first frame was skipped
597 file.write('set cframes [expr $nframes - 1]\n')
598 file.write('puts $nframes_file $cframes\n')
599 # Remove duplicated indices
600 # Use the -integer argument to do the correct sorting
601 # Otherwise numbers are sorted as strings so you have 1, 10, 100, 2, etc.
602 file.write('set unique_interface1_atom_indices [lsort -integer -unique $accumulated_interface1_atom_indices]\n')
603 file.write('set unique_interface2_atom_indices [lsort -integer -unique $accumulated_interface2_atom_indices]\n')
604 # Write indices to files
605 file.write(f'set indices_file [open {interface_selection_1_filename} w]\n')
606 file.write('puts $indices_file $unique_interface1_atom_indices\n')
607 file.write(f'set indices_file [open {interface_selection_2_filename} w]\n')
608 file.write('puts $indices_file $unique_interface2_atom_indices\n')
609 file.write('exit\n')
611 # Run VMD
612 logs = run([
613 "vmd",
614 input_structure_filepath,
615 input_trajectory_filepath,
616 "-e",
617 commands_filename,
618 "-dispdev",
619 "none"
620 ], stdout=PIPE, stderr=PIPE).stdout.decode()
622 # If any of the output files do not exist at this point then it means something went wrong with vmd
623 expected_output_files = [
624 selection_1_filename,
625 selection_2_filename,
626 interface_selection_1_filename,
627 interface_selection_2_filename,
628 interacting_frames_filename,
629 total_frames_filename
630 ]
631 for output_file in expected_output_files:
632 if not os.path.exists(output_file):
633 print(logs)
634 raise SystemExit('Something went wrong with VMD')
636 # Set a function to read the VMD output and parse the atom indices string to an array of integers
637 def process_vmd_output (output_filename : str) -> list[int]:
638 with open(output_filename, 'r') as file:
639 raw_atom_indices = file.read()
640 return [ int(i) for i in raw_atom_indices.split() ]
642 # Read the VMD output
643 selection_1_atom_indices = process_vmd_output(selection_1_filename)
644 selection_2_atom_indices = process_vmd_output(selection_2_filename)
645 selection_1_interface_atom_indices = process_vmd_output(interface_selection_1_filename)
646 selection_2_interface_atom_indices = process_vmd_output(interface_selection_2_filename)
647 interacting_frames = process_vmd_output(interacting_frames_filename)[0]
648 total_frames = process_vmd_output(total_frames_filename)[0]
650 # Remove trash files
651 trash_files = [ commands_filename ] + expected_output_files
652 for trash_file in trash_files:
653 os.remove(trash_file)
655 # Return the results
656 return {
657 'selection_1_atom_indices': selection_1_atom_indices,
658 'selection_2_atom_indices': selection_2_atom_indices,
659 'selection_1_interface_atom_indices': selection_1_interface_atom_indices,
660 'selection_2_interface_atom_indices': selection_2_interface_atom_indices,
661 'interacting_frames': interacting_frames,
662 'total_frames': total_frames
663 }