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

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/ 

4 

5import os 

6from os.path import exists 

7 

8from subprocess import run, PIPE, STDOUT 

9 

10from model_workflow.utils.file import File 

11from model_workflow.utils.type_hints import * 

12from model_workflow.utils.auxiliar import warn 

13 

14# Set characters to be escaped since they have a meaning in TCL 

15TCL_RESERVED_CHARACTERS = ['"','[',']'] 

16 

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 

23 

24# Set the script filename with all commands to be passed to vmd 

25commands_filename = '.commands.vmd' 

26 

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'} 

30 

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} 

47 

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): 

53 

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] 

58 

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') 

69 

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() 

80 

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') 

86 

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] 

100 

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): 

117 

118 # If no atom selection is provided then all atoms are selected 

119 if not atom_selection: 

120 atom_selection = 'all' 

121 

122 # Escape TCL meaningful characters 

123 escaped_atom_selection = escape_tcl_selection(atom_selection) 

124 

125 # If no chain letter is provided then the flag 'fragment' is used 

126 if not chain_letter: 

127 chain_letter = 'fragment' 

128 

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 

132 

133 # Check the file exists 

134 if not exists(input_pdb_filename): 

135 raise SystemExit('ERROR: The file does not exist') 

136 

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') 

163 

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() 

173 

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') 

178 

179 # Remove the vmd commands file 

180 os.remove(commands_filename) 

181 

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 ): 

190 

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') 

193 

194 # Get the format to export coordinates 

195 output_trajectory_file = File(output_trajectory_filename) 

196 output_trajectory_format = output_trajectory_file.format 

197 

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' 

201 

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') 

211 

212 inputs = [ input_structure_filename, *input_trajectory_filenames ] if input_structure_filename else input_trajectory_filenames 

213 

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 

223 

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') 

228 

229 # Remove the vmd commands file 

230 os.remove(commands_filename) 

231 

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] 

255 

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]: 

258 

259 # Escape TCL meaningful characters 

260 escaped_selection = escape_tcl_selection(selection) 

261 

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') 

274 

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() 

284 

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') 

289 

290 # Read the VMD output 

291 with open(atom_indices_filename, 'r') as file: 

292 raw_atom_indices = file.read() 

293 

294 # Parse the atom indices string to an array of integers 

295 atom_indices = [ int(i) for i in raw_atom_indices.split() ] 

296 

297 # Remove trahs files 

298 trash_files = [ commands_filename, atom_indices_filename ] 

299 for trash_file in trash_files: 

300 os.remove(trash_file) 

301 

302 return atom_indices 

303 

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] ]: 

307 

308 # Parse the selection to vmd 

309 vmd_selection = 'all' 

310 if selection: 

311 vmd_selection = selection.to_vmd() 

312 

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') 

324 

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() 

334 

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') 

339 

340 # Read the VMD output 

341 with open(output_bonds_file, 'r') as file: 

342 raw_bonds = file.read() 

343 

344 # Remove vmd files since they are no longer usefull 

345 for f in [ commands_filename, output_bonds_file ]: 

346 os.remove(f) 

347 

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', '') + ' ' 

351 

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 

385 

386 return bonds_per_atom 

387 

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] ]: 

395 

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() 

399 

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') 

425 

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() 

435 

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') 

440 

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() 

448 

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) 

452 

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', '') + ' ' 

456 

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() ] 

460 

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 

494 

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') 

498 

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) 

507 

508 return crossed_bonds 

509 

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]: 

520 

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})') 

524 

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' 

533 

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') 

605 

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() 

616 

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') 

630 

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() ] 

636 

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] 

644 

645 # Remove trash files 

646 trash_files = [ commands_filename ] + expected_output_files 

647 for trash_file in trash_files: 

648 os.remove(trash_file) 

649 

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 }