1 #!/bin/env python
2 #
3 # File: VinaPerformDocking.py
4 # Author: Manish Sud <msud@san.rr.com>
5 #
6 # Acknowledgments: Diogo Santos-Martins and Stefano Forli
7 #
8 # Copyright (C) 2026 Manish Sud. All rights reserved.
9 #
10 # The functionality available in this script is implemented using AutoDockVina
11 # and Meeko, open source software packages for docking, and RDKit, an open
12 # source toolkit for cheminformatics developed by Greg Landrum.
13 #
14 # This file is part of MayaChemTools.
15 #
16 # MayaChemTools is free software; you can redistribute it and/or modify it under
17 # the terms of the GNU Lesser General Public License as published by the Free
18 # Software Foundation; either version 3 of the License, or (at your option) any
19 # later version.
20 #
21 # MayaChemTools is distributed in the hope that it will be useful, but without
22 # any warranty; without even the implied warranty of merchantability of fitness
23 # for a particular purpose. See the GNU Lesser General Public License for more
24 # details.
25 #
26 # You should have received a copy of the GNU Lesser General Public License
27 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
28 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
29 # Boston, MA, 02111-1307, USA.
30 #
31
32 from __future__ import print_function
33
34 import os
35 import sys
36 import time
37 import re
38 import tempfile
39 import json
40 import glob
41 import multiprocessing as mp
42
43 # AutoDock Vina imports...
44 try:
45 from vina import Vina
46 from vina import __version__ as vinaVersion
47 except ImportError as ErrMsg:
48 sys.stderr.write("\nFailed to import AutoDock Vina module/package: %s\n" % ErrMsg)
49 sys.stderr.write("Check/update your Vina environment and try again.\n\n")
50 sys.exit(1)
51
52 # AutoDock Meeko imports...
53 try:
54 from meeko import __version__ as meekoVersion
55 from meeko import MoleculePreparation
56 from meeko import PDBQTMolecule
57 from meeko import RDKitMolCreate
58 from meeko import PDBQTWriterLegacy
59 except ImportError as ErrMsg:
60 sys.stderr.write("\nFailed to import AutoDock Meeko module/package: %s\n" % ErrMsg)
61 sys.stderr.write("Check/update your Meeko environment and try again.\n\n")
62 sys.exit(1)
63
64 # RDKit imports...
65 try:
66 from rdkit import rdBase
67 from rdkit import Chem
68 from rdkit.Chem import rdMolTransforms
69 except ImportError as ErrMsg:
70 sys.stderr.write("\nFailed to import RDKit module/package: %s\n" % ErrMsg)
71 sys.stderr.write("Check/update your RDKit environment and try again.\n\n")
72 sys.exit(1)
73
74 # MayaChemTools imports...
75 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
76 try:
77 from docopt import docopt
78 import MiscUtil
79 import RDKitUtil
80 except ImportError as ErrMsg:
81 sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
82 sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
83 sys.exit(1)
84
85 ScriptName = os.path.basename(sys.argv[0])
86 Options = {}
87 OptionsInfo = {}
88
89
90 def main():
91 """Start execution of the script."""
92
93 MiscUtil.PrintInfo(
94 "\n%s (Vina v%s; Meeko v%s; RDKit v%s; MayaChemTools v%s; %s): Starting...\n"
95 % (
96 ScriptName,
97 vinaVersion,
98 meekoVersion,
99 rdBase.rdkitVersion,
100 MiscUtil.GetMayaChemToolsVersion(),
101 time.asctime(),
102 )
103 )
104
105 (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
106
107 # Retrieve command line arguments and options...
108 RetrieveOptions()
109
110 # Process and validate command line arguments and options...
111 ProcessOptions()
112
113 # Perform actions required by the script...
114 PerformDocking()
115
116 MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
117 MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
118
119
120 def PerformDocking():
121 """Perform docking."""
122
123 # Setup a molecule reader...
124 MiscUtil.PrintInfo("\nProcessing file %s..." % OptionsInfo["Infile"])
125 Mols = RDKitUtil.ReadMolecules(OptionsInfo["Infile"], **OptionsInfo["InfileParams"])
126
127 # Set up molecule writers...
128 Writer, WriterFlexRes = SetupWriters()
129 if WriterFlexRes is None:
130 MiscUtil.PrintInfo("Generating file %s..." % OptionsInfo["Outfile"])
131 else:
132 MiscUtil.PrintInfo("Generating files %s and %s..." % (OptionsInfo["Outfile"], OptionsInfo["OutfileFlexRes"]))
133
134 MolCount, ValidMolCount, DockingFailedCount = ProcessMolecules(Mols, Writer, WriterFlexRes)
135
136 CloseWriters(Writer, WriterFlexRes)
137
138 MiscUtil.PrintInfo("\nTotal number of molecules: %d" % MolCount)
139 MiscUtil.PrintInfo("Number of valid molecules: %d" % ValidMolCount)
140 MiscUtil.PrintInfo("Number of molecules failed during docking: %d" % DockingFailedCount)
141 MiscUtil.PrintInfo("Number of ignored molecules: %d" % (MolCount - ValidMolCount + DockingFailedCount))
142
143
144 def ProcessMolecules(Mols, Writer, WriterFlexRes):
145 """Process and dock molecules."""
146
147 if OptionsInfo["MPMode"]:
148 return ProcessMoleculesUsingMultipleProcesses(Mols, Writer, WriterFlexRes)
149 else:
150 return ProcessMoleculesUsingSingleProcess(Mols, Writer, WriterFlexRes)
151
152
153 def ProcessMoleculesUsingSingleProcess(Mols, Writer, WriterFlexRes):
154 """Process and dock molecules using a single process."""
155
156 VinaHandle = InitializeVina(OptionsInfo["QuietMode"])
157 MolPrepHandle = InitializeMeekoMolPrepration(OptionsInfo["QuietMode"])
158 if not OptionsInfo["QuietMode"]:
159 MiscUtil.PrintInfo("\nVina configuration:\n%s" % VinaHandle)
160
161 MiscUtil.PrintInfo("\nDocking molecules (Mode: %s)..." % OptionsInfo["Mode"])
162
163 (MolCount, ValidMolCount, DockingFailedCount) = [0] * 3
164 for Mol in Mols:
165 MolCount += 1
166
167 if Mol is None:
168 continue
169
170 if not CheckAndValidateMolecule(Mol, MolCount):
171 continue
172
173 ValidMolCount += 1
174 CalcStatus, PoseMol, PoseMolEnergies, PoseMolFlexRes = DockMolecule(VinaHandle, MolPrepHandle, Mol, MolCount)
175
176 if not CalcStatus:
177 if not OptionsInfo["QuietMode"]:
178 MiscUtil.PrintWarning("Failed to dock molecule %s" % RDKitUtil.GetMolName(Mol, MolCount))
179
180 DockingFailedCount += 1
181 continue
182
183 WriteMolPoses(Writer, Mol, MolCount, PoseMol, PoseMolEnergies, WriterFlexRes, PoseMolFlexRes)
184
185 return (MolCount, ValidMolCount, DockingFailedCount)
186
187
188 def ProcessMoleculesUsingMultipleProcesses(Mols, Writer, WriterFlexRes):
189 """Process and calculate energy of molecules using multiprocessing."""
190
191 MiscUtil.PrintInfo("\nDocking molecules using multiprocessing...")
192 MiscUtil.PrintInfo("\nDocking molecules (Mode: %s)..." % OptionsInfo["Mode"])
193
194 MPParams = OptionsInfo["MPParams"]
195
196 # Setup data for initializing a worker process...
197 InitializeWorkerProcessArgs = (
198 MiscUtil.ObjectToBase64EncodedString(Options),
199 MiscUtil.ObjectToBase64EncodedString(OptionsInfo),
200 )
201
202 # Setup a encoded mols data iterable for a worker process...
203 WorkerProcessDataIterable = RDKitUtil.GenerateBase64EncodedMolStrings(Mols)
204
205 # Setup process pool along with data initialization for each process...
206 if not OptionsInfo["QuietMode"]:
207 MiscUtil.PrintInfo(
208 "\nConfiguring multiprocessing using %s method..."
209 % ("mp.Pool.imap()" if re.match("^Lazy$", MPParams["InputDataMode"], re.I) else "mp.Pool.map()")
210 )
211 MiscUtil.PrintInfo(
212 "NumProcesses: %s; InputDataMode: %s; ChunkSize: %s\n"
213 % (
214 MPParams["NumProcesses"],
215 MPParams["InputDataMode"],
216 ("automatic" if MPParams["ChunkSize"] is None else MPParams["ChunkSize"]),
217 )
218 )
219
220 ProcessPool = mp.Pool(MPParams["NumProcesses"], InitializeWorkerProcess, InitializeWorkerProcessArgs)
221
222 # Start processing...
223 if re.match("^Lazy$", MPParams["InputDataMode"], re.I):
224 Results = ProcessPool.imap(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
225 elif re.match("^InMemory$", MPParams["InputDataMode"], re.I):
226 Results = ProcessPool.map(WorkerProcess, WorkerProcessDataIterable, MPParams["ChunkSize"])
227 else:
228 MiscUtil.PrintError(
229 'The value, %s, specified for "--inputDataMode" is not supported.' % (MPParams["InputDataMode"])
230 )
231
232 (MolCount, ValidMolCount, DockingFailedCount) = [0] * 3
233 for Result in Results:
234 MolCount += 1
235
236 MolIndex, EncodedMol, CalcStatus, EncodedPoseMol, PoseMolEnergies, EncodedPoseMolFlexRes = Result
237
238 if EncodedMol is None:
239 continue
240
241 ValidMolCount += 1
242
243 if not CalcStatus:
244 if not OptionsInfo["QuietMode"]:
245 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
246 MolName = RDKitUtil.GetMolName(Mol, MolCount)
247 MiscUtil.PrintWarning("Failed to dock molecule %s" % MolName)
248
249 DockingFailedCount += 1
250 continue
251
252 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
253 PoseMol = None if EncodedPoseMol is None else RDKitUtil.MolFromBase64EncodedMolString(EncodedPoseMol)
254 PoseMolFlexRes = (
255 None if EncodedPoseMolFlexRes is None else RDKitUtil.MolFromBase64EncodedMolString(EncodedPoseMolFlexRes)
256 )
257
258 WriteMolPoses(Writer, Mol, MolCount, PoseMol, PoseMolEnergies, WriterFlexRes, PoseMolFlexRes)
259
260 return (MolCount, ValidMolCount, DockingFailedCount)
261
262
263 def InitializeWorkerProcess(*EncodedArgs):
264 """Initialize data for a worker process."""
265
266 global Options, OptionsInfo
267
268 if not OptionsInfo["QuietMode"]:
269 MiscUtil.PrintInfo("Starting process (PID: %s)..." % os.getpid())
270
271 # Decode Options and OptionInfo...
272 Options = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[0])
273 OptionsInfo = MiscUtil.ObjectFromBase64EncodedString(EncodedArgs[1])
274
275 # Initialize in a quiet mode...
276 OptionsInfo["VinaHandle"] = InitializeVina(True)
277 OptionsInfo["MolPrepHandle"] = InitializeMeekoMolPrepration(True)
278
279
280 def WorkerProcess(EncodedMolInfo):
281 """Process data for a worker process."""
282
283 MolIndex, EncodedMol = EncodedMolInfo
284
285 CalcStatus, PoseMol, PoseMolEnergies, PoseMolFlexRes = [False, None, None, None]
286
287 if EncodedMol is None:
288 return [MolIndex, None, CalcStatus, PoseMol, PoseMolEnergies, PoseMolFlexRes]
289
290 Mol = RDKitUtil.MolFromBase64EncodedMolString(EncodedMol)
291 MolCount = MolIndex + 1
292
293 if not CheckAndValidateMolecule(Mol, MolCount):
294 return [MolIndex, None, CalcStatus, PoseMol, PoseMolEnergies, PoseMolFlexRes]
295
296 CalcStatus, PoseMol, PoseMolEnergies, PoseMolFlexRes = DockMolecule(
297 OptionsInfo["VinaHandle"], OptionsInfo["MolPrepHandle"], Mol, MolCount
298 )
299
300 EncodedPoseMol = (
301 None
302 if PoseMol is None
303 else RDKitUtil.MolToBase64EncodedMolString(
304 PoseMol, PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps
305 )
306 )
307
308 EncodedPoseMolFlexRes = (
309 None
310 if PoseMolFlexRes is None
311 else RDKitUtil.MolToBase64EncodedMolString(
312 PoseMolFlexRes,
313 PropertyPickleFlags=Chem.PropertyPickleOptions.MolProps | Chem.PropertyPickleOptions.PrivateProps,
314 )
315 )
316
317 return [MolIndex, EncodedMol, CalcStatus, EncodedPoseMol, PoseMolEnergies, EncodedPoseMolFlexRes]
318
319
320 def DockMolecule(VinaHandle, MolPrepHandle, Mol, MolNum=None):
321 """Dock molecule."""
322
323 Status, PoseMol, PoseMolEnergies, PoseMolFlexRes = [False, None, None, None]
324
325 if OptionsInfo["VinaVerbosity"] != 0:
326 MiscUtil.PrintInfo("\nProcessing molecule %s..." % RDKitUtil.GetMolName(Mol, MolNum))
327
328 PDBQTMolStr = PrepareMolecule(MolPrepHandle, Mol)
329 if PDBQTMolStr is None:
330 return (Status, PoseMol, PoseMolEnergies, PoseMolFlexRes)
331 VinaHandle.set_ligand_from_string(PDBQTMolStr)
332
333 if OptionsInfo["ScoreOnlyMode"]:
334 return ProcessMoleculeForScoreOnlyMode(VinaHandle, Mol)
335 elif OptionsInfo["LocalOptimizationOnlyMode"]:
336 return ProcessMoleculeForLocalOptimizationOnlyMode(VinaHandle, Mol)
337 elif OptionsInfo["DockMode"]:
338 return ProcessMoleculeForDockMode(VinaHandle, Mol)
339 else:
340 MiscUtil.PrintError(
341 'The value specified, %s, for option "-m, --mode" is not valid. Supported values: Dock LocalOptimizationOnly ScoreOnly'
342 % OptionsInfo["Mode"]
343 )
344
345 return (Status, PoseMol, PoseMolEnergies, PoseMolFlexRes)
346
347
348 def ProcessMoleculeForScoreOnlyMode(VinaHandle, Mol):
349 """Score molecule without performing any docking."""
350
351 Status, PoseMol, Energies, PoseMolFlexRes = [False, None, None, None]
352
353 # Score molecule...
354 try:
355 Energies = VinaHandle.score()
356 except Exception as ErrMsg:
357 if not OptionsInfo["QuietMode"]:
358 MiscUtil.PrintWarning("Failed to score molecule:\n%s\n" % ErrMsg)
359 return (False, None, None, None)
360
361 if len(Energies) == 0:
362 return (False, None, None, None)
363
364 Status = True
365 PoseMol = Mol
366
367 return (Status, PoseMol, Energies, PoseMolFlexRes)
368
369
370 def ProcessMoleculeForLocalOptimizationOnlyMode(VinaHandle, Mol):
371 """Score molecule after a local optimization wthout any docking."""
372
373 Status, PoseMol, Energies, PoseMolFlexRes = [False, None, None, None]
374
375 # Optimize and score molecule...
376 try:
377 Energies = VinaHandle.optimize()
378 except Exception as ErrMsg:
379 if not OptionsInfo["QuietMode"]:
380 MiscUtil.PrintWarning("Failed to perform local optimization:\n%s\n" % ErrMsg)
381 return (False, None, None, None)
382
383 if len(Energies) == 0:
384 return (False, None, None, None)
385
386 # Write optimize pose to a temporray file and retrieve the PDBQT string for the pose...
387 (_, TmpFile) = tempfile.mkstemp(suffix=".pdbqt", prefix="VinaOptimize_", text=True)
388 VinaHandle.write_pose(TmpFile, overwrite=True)
389
390 TmpFH = open(TmpFile, "r")
391 PosePDBQTOutputStr = TmpFH.read()
392 TmpFH.close()
393
394 os.remove(TmpFile)
395
396 # Setup a mol containing optimized pose...
397 try:
398 PosePDBQTOutputMol = PDBQTMolecule(PosePDBQTOutputStr)
399 PoseMols = RDKitMolCreate.from_pdbqt_mol(PosePDBQTOutputMol)
400 except Exception as ErrMsg:
401 if not OptionsInfo["QuietMode"]:
402 MiscUtil.PrintWarning("Failed to retrieve optimize pose:\n%s\n" % ErrMsg)
403 return (False, None, None, None)
404
405 if len(PoseMols) == 0:
406 return (False, None, None, None)
407
408 Status = True
409 PoseMol = PoseMols[0]
410
411 return (Status, PoseMol, Energies, PoseMolFlexRes)
412
413
414 def ProcessMoleculeForDockMode(VinaHandle, Mol):
415 """Dock and score molecule."""
416
417 Status, PoseMol, Energies, PoseMolFlexRes = [False, None, None, None]
418
419 # Dock molecule...
420 try:
421 VinaHandle.dock(
422 exhaustiveness=OptionsInfo["Exhaustiveness"],
423 n_poses=OptionsInfo["NumPoses"],
424 min_rmsd=OptionsInfo["MinRMSD"],
425 max_evals=OptionsInfo["MaxEvaluations"],
426 )
427 except Exception as ErrMsg:
428 if not OptionsInfo["QuietMode"]:
429 MiscUtil.PrintWarning("Failed to dock molecule:\n%s\n" % ErrMsg)
430 return (False, None, None, None)
431
432 # Retrieve poses...
433 try:
434 PosePDBQTOutputStr = VinaHandle.poses(
435 n_poses=OptionsInfo["NumPoses"], energy_range=OptionsInfo["EnergyRange"], coordinates_only=False
436 )
437 except Exception as ErrMsg:
438 if not OptionsInfo["QuietMode"]:
439 MiscUtil.PrintWarning("Failed to retrieve docked poses:\n%s\n" % ErrMsg)
440 return (False, None, None, None)
441
442 PoseMol, PoseMolFlexRes = SetupPoseMols(PosePDBQTOutputStr)
443 if PoseMol is None:
444 return (False, None, None, None)
445
446 # Retrieve energies...
447 try:
448 Energies = VinaHandle.energies(n_poses=OptionsInfo["NumPoses"], energy_range=OptionsInfo["EnergyRange"])
449 except Exception as ErrMsg:
450 if not OptionsInfo["QuietMode"]:
451 MiscUtil.PrintWarning("Failed to retrieve energies for docked poses:\n%s\n" % ErrMsg)
452 return (False, None, None, None)
453
454 if len(Energies) == 0:
455 return (False, None, None, None)
456
457 Status = True
458
459 return (Status, PoseMol, Energies, PoseMolFlexRes)
460
461
462 def SetupPoseMols(PosePDBQTOutputStr):
463 """Process PDBQT Vina poses string to setup pose mols for the docked
464 molecule and any flexible residues."""
465
466 PoseMol, PoseMolFlexRes = [None, None]
467
468 # Setup a mol containing poses as conformers...
469 try:
470 PosePDBQTOutputMol = PDBQTMolecule(PosePDBQTOutputStr)
471 PoseMols = RDKitMolCreate.from_pdbqt_mol(PosePDBQTOutputMol)
472 except Exception as ErrMsg:
473 if not OptionsInfo["QuietMode"]:
474 MiscUtil.PrintWarning("Failed to retrieve docked poses:\n%s\n" % ErrMsg)
475 return (None, None)
476
477 if len(PoseMols) == 0:
478 return (None, None)
479
480 # First mol is the pose for docked molecule...
481 PoseMol = PoseMols.pop(0)
482
483 # Collect pose mols foe flexible side chain residues...
484 PoseMolsFlexRes = []
485 for Mol in PoseMols:
486 if not Mol.HasProp("meeko"):
487 continue
488 MeekoPropMap = json.loads(Mol.GetProp("meeko"))
489 if type(MeekoPropMap) is dict:
490 if "is_sidechain" in MeekoPropMap:
491 if MeekoPropMap["is_sidechain"]:
492 PoseMolsFlexRes.append(Mol)
493
494 # Combine pose mols for flexible side chain residues into a single pose mol...
495 if len(PoseMolsFlexRes):
496 for Mol in PoseMolsFlexRes:
497 if PoseMolFlexRes is None:
498 PoseMolFlexRes = Mol
499 else:
500 PoseMolFlexRes = Chem.CombineMols(PoseMolFlexRes, Mol)
501
502 if PoseMolFlexRes is not None:
503 PoseMolConfCount = PoseMol.GetNumConformers()
504 PoseMolFlexResConfCount = PoseMolFlexRes.GetNumConformers()
505 if PoseMolConfCount != PoseMolFlexResConfCount:
506 if not OptionsInfo["QuietMode"]:
507 MiscUtil.PrintWarning(
508 "The number of poses, %s, for flexible residues doesn't match number of poses, %s, for docked molecule...\n"
509 % (PoseMolFlexResConfCount, PoseMolConfCount)
510 )
511
512 return (PoseMol, PoseMolFlexRes)
513
514
515 def PrepareMolecule(MolPrepHandle, Mol):
516 """Prepare molecule for docking."""
517
518 try:
519 PreppedMols = MolPrepHandle.prepare(Mol)
520 except Exception as ErrMsg:
521 if not OptionsInfo["QuietMode"]:
522 MiscUtil.PrintWarning("Failed to prepare molecule for docking:\n%s\n" % ErrMsg)
523 return None
524
525 if len(PreppedMols) == 0:
526 return None
527
528 PreppedMol = PreppedMols[0]
529
530 # Setup PDBQT mole string...
531 try:
532 PDBQTMolStr, Status, ErrMsg = PDBQTWriterLegacy.write_string(PreppedMol)
533 except Exception as ExceptionErrMsg:
534 if not OptionsInfo["QuietMode"]:
535 MiscUtil.PrintWarning("Failed to prepare molecule for docking:\n%s\n" % ExceptionErrMsg)
536 return None
537
538 if not Status:
539 if not OptionsInfo["QuietMode"]:
540 MiscUtil.PrintWarning("Failed to prepare molecule for docking:\n%s\n" % ErrMsg)
541 return None
542
543 if MiscUtil.IsEmpty(PDBQTMolStr):
544 return None
545
546 return PDBQTMolStr
547
548
549 def InitializeVina(Quiet=False):
550 """Initialize AutoDock Vina."""
551
552 if not Quiet:
553 MiscUtil.PrintInfo("\nInitializing Vina...")
554
555 VinaHandle = Vina(
556 sf_name=OptionsInfo["Forcefield"],
557 cpu=OptionsInfo["NumThreads"],
558 seed=OptionsInfo["RandomSeed"],
559 no_refine=OptionsInfo["SkipRefinement"],
560 verbosity=OptionsInfo["VinaVerbosity"],
561 )
562
563 SetupReceptor(VinaHandle, Quiet)
564 SetupForcefieldWeights(VinaHandle, Quiet)
565 SetupReceptorMaps(VinaHandle, Quiet)
566
567 return VinaHandle
568
569
570 def SetupReceptor(VinaHandle, Quiet=False):
571 """Setup receptor"""
572
573 if OptionsInfo["UseReceptorFile"] and OptionsInfo["UseReceptorFlexFile"]:
574 VinaHandle.set_receptor(
575 rigid_pdbqt_filename=OptionsInfo["ReceptorFile"], flex_pdbqt_filename=OptionsInfo["ReceptorFlexFile"]
576 )
577 elif OptionsInfo["UseReceptorFile"]:
578 VinaHandle.set_receptor(rigid_pdbqt_filename=OptionsInfo["ReceptorFile"], flex_pdbqt_filename=None)
579 elif OptionsInfo["UseReceptorFlexFile"]:
580 VinaHandle.set_receptor(rigid_pdbqt_filename=None, flex_pdbqt_filename=OptionsInfo["ReceptorFlexFile"])
581
582
583 def SetupForcefieldWeights(VinaHandle, Quiet=False):
584 """Setup forcefield weights."""
585
586 Weights = OptionsInfo["ForcefieldWeightParams"]
587 Forcefield = OptionsInfo["Forcefield"]
588 if not OptionsInfo["ForcefieldWeightParamsSpecified"]:
589 if not Quiet:
590 MiscUtil.PrintInfo("\nUsing default forcefield weights for %s..." % Forcefield)
591 return
592
593 if not Quiet:
594 MiscUtil.PrintInfo("\nSetting specified forcefield weights for %s..." % Forcefield)
595
596 if OptionsInfo["UseAD4Forcefield"]:
597 VinaHandle.set_weights(
598 [
599 Weights["AD4Vdw"],
600 Weights["AD4HydrogenBond"],
601 Weights["AD4Electrostatic"],
602 Weights["AD4Desolvation"],
603 Weights["AD4GlueLinearAttraction"],
604 Weights["AD4Rot"],
605 ]
606 )
607 elif OptionsInfo["UseVinaForcefield"]:
608 VinaHandle.set_weights(
609 [
610 Weights["VinaGaussian1"],
611 Weights["VinaGaussian2"],
612 Weights["VinaRepulsion"],
613 Weights["VinaHydrophobic"],
614 Weights["VinaHydrogenBond"],
615 Weights["VinaGlueLinearAttraction"],
616 Weights["VinaRot"],
617 ]
618 )
619 elif OptionsInfo["UseVinardoForcefield"]:
620 VinaHandle.set_weights(
621 [
622 Weights["VinardoGaussian1"],
623 Weights["VinardoRepulsion"],
624 Weights["VinardoHydrophobic"],
625 Weights["VinardoHydrogenBond"],
626 Weights["VinardoGlueLinearAttraction"],
627 Weights["VinardoRot"],
628 ]
629 )
630 else:
631 MiscUtil.PrintError(
632 'The value specified, %s, for option "-f, --forcefield" is not valid. Supported values: AD4 Vina Vinardo'
633 % Forcefield
634 )
635
636
637 def SetupReceptorMaps(VinaHandle, Quiet=False):
638 """Setup receptor maps."""
639
640 if not Quiet:
641 MiscUtil.PrintInfo("\nSetting up receptor and maps for %s forcefield..." % OptionsInfo["Forcefield"])
642
643 if OptionsInfo["UseAD4Forcefield"]:
644 # Load maps for rigid portion of the receptor...
645 VinaHandle.load_maps(OptionsInfo["ReceptorMapsPrefix"])
646 elif OptionsInfo["UseVinaForcefield"] or OptionsInfo["UseVinardoForcefield"]:
647 # Load or compute maps for rigid portion of the receptor...
648 if OptionsInfo["UseReceptorMapsPrefix"]:
649 VinaHandle.load_maps(OptionsInfo["ReceptorMapsPrefix"])
650 else:
651 VinaHandle.compute_vina_maps(
652 center=OptionsInfo["GridCenterList"],
653 box_size=OptionsInfo["GridSizeList"],
654 spacing=OptionsInfo["GridSpacing"],
655 force_even_voxels=False,
656 )
657 else:
658 MiscUtil.PrintError(
659 'The value specified, %s, for option "-f, --forcefield" is not valid. Supported values: AD4 Vina Vinardo'
660 % OptionsInfo["Forcefield"]
661 )
662
663
664 def InitializeMeekoMolPrepration(Quiet=False):
665 """Initialize meeko molecule prepration."""
666
667 if OptionsInfo["MergeHydrogens"]:
668 MolPrep = MoleculePreparation(merge_these_atom_types=("H",))
669 else:
670 MolPrep = MoleculePreparation(merge_these_atom_types="")
671
672 return MolPrep
673
674
675 def CheckAndValidateMolecule(Mol, MolCount=None):
676 """Validate molecule for docking."""
677
678 MolName = RDKitUtil.GetMolName(Mol, MolCount)
679 if RDKitUtil.IsMolEmpty(Mol):
680 if not OptionsInfo["QuietMode"]:
681 MiscUtil.PrintWarning("Ignoring empty molecule: %s\n" % MolName)
682 return False
683
684 if not OptionsInfo["ValidateMolecules"]:
685 return True
686
687 # Check for 3D flag...
688 if not Mol.GetConformer().Is3D():
689 if not OptionsInfo["QuietMode"]:
690 MiscUtil.PrintWarning("3D tag is not set for molecule: %s\n" % MolName)
691
692 # Check for missing hydrogens...
693 if RDKitUtil.AreHydrogensMissingInMolecule(Mol):
694 if not OptionsInfo["QuietMode"]:
695 MiscUtil.PrintWarning("Missing hydrogens in molecule: %s\n" % MolName)
696
697 return True
698
699
700 def WriteMolPoses(Writer, Mol, MolNum, PoseMol, Energies, WriterFlexRes=None, PoseMolFlexRes=None):
701 """Write molecule."""
702
703 if OptionsInfo["ScoreOnlyMode"]:
704 return WriteMolPoseForScoreOnlyMode(Writer, Mol, MolNum, PoseMol, Energies)
705 elif OptionsInfo["LocalOptimizationOnlyMode"]:
706 return WriteMolPoseForLocalOptimizationOnlyMode(Writer, Mol, MolNum, PoseMol, Energies)
707 elif OptionsInfo["DockMode"]:
708 return WriteMolPosesForDockMode(Writer, Mol, MolNum, PoseMol, Energies, WriterFlexRes, PoseMolFlexRes)
709 else:
710 MiscUtil.PrintError(
711 'The value specified, %s, for option "-m, --mode" is not valid. Supported values: Dock LocalOptimizationOnly ScoreOnly'
712 % OptionsInfo["Mode"]
713 )
714
715
716 def WriteMolPoseForScoreOnlyMode(Writer, Mol, MolNum, PoseMol, Energies):
717 """Write out molecule and associated information for score only mode."""
718
719 ClearMeekoMolProperties(PoseMol)
720
721 MolName = RDKitUtil.GetMolName(Mol, MolNum)
722 PoseMol.SetProp("_Name", MolName)
723
724 SetupEnergyProperties(PoseMol, Energies)
725 Writer.write(PoseMol)
726
727 return
728
729
730 def WriteMolPoseForLocalOptimizationOnlyMode(Writer, Mol, MolNum, PoseMol, Energies):
731 """Write out molecule and associated information for score only mode."""
732
733 ClearMeekoMolProperties(PoseMol)
734
735 MolName = RDKitUtil.GetMolName(Mol, MolNum)
736 PoseMol.SetProp("_Name", MolName)
737
738 SetupEnergyProperties(PoseMol, Energies)
739 Writer.write(PoseMol)
740
741
742 def WriteMolPosesForDockMode(Writer, Mol, MolNum, PoseMol, Energies, WriterFlexRes=None, PoseMolFlexRes=None):
743 """Write out molecule and associated information for dock mode."""
744
745 MolName = RDKitUtil.GetMolName(Mol, MolNum)
746
747 # Write out poses for docked molecule...
748 ClearMeekoMolProperties(PoseMol)
749 for PoseMolConfIndex, PoseMolConf in enumerate(PoseMol.GetConformers()):
750 PoseMolName = "%s_Pose%s" % (MolName, (PoseMolConfIndex + 1))
751 PoseMol.SetProp("_Name", PoseMolName)
752
753 SetupEnergyProperties(PoseMol, Energies[PoseMolConfIndex])
754
755 Writer.write(PoseMol, confId=PoseMolConf.GetId())
756
757 # Write out poses for flexible reside side chains...
758 if WriterFlexRes is None or PoseMolFlexRes is None:
759 return
760
761 ClearMeekoMolProperties(PoseMolFlexRes)
762 for PoseMolFlexResConfIndex, PoseMolFlexResConf in enumerate(PoseMolFlexRes.GetConformers()):
763 PoseMolFlexResName = "%s_Flex_Receptor_Pose%s" % (MolName, (PoseMolFlexResConfIndex + 1))
764 PoseMolFlexRes.SetProp("_Name", PoseMolFlexResName)
765
766 WriterFlexRes.write(PoseMolFlexRes, confId=PoseMolFlexResConf.GetId())
767
768
769 def ClearMeekoMolProperties(Mol):
770 """Clear Meeko molecule properties."""
771
772 for PropName in ["meeko"]:
773 if Mol.HasProp(PropName):
774 Mol.ClearProp(PropName)
775
776
777 def SetupEnergyProperties(Mol, Energies):
778 """Setup energy properties."""
779
780 Precision = OptionsInfo["Precision"]
781
782 for Index, EnergyLabel in enumerate(OptionsInfo["EnergyLabelsList"]):
783 EnergyValueIndex = OptionsInfo["EnergyValueIndicesList"][Index]
784 EnergyValue = "%.*f" % (Precision, Energies[EnergyValueIndex])
785 Mol.SetProp(EnergyLabel, EnergyValue)
786
787
788 def SetupWriters():
789 """Setup writers for output files."""
790
791 Writer, WriterFlexRes = [None, None]
792
793 Writer = RDKitUtil.MoleculesWriter(OptionsInfo["Outfile"], **OptionsInfo["OutfileParams"])
794 if Writer is None:
795 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["Outfile"])
796
797 if OptionsInfo["DockMode"] and OptionsInfo["UseReceptorFlexFile"]:
798 WriterFlexRes = RDKitUtil.MoleculesWriter(OptionsInfo["OutfileFlexRes"], **OptionsInfo["OutfileParams"])
799 if WriterFlexRes is None:
800 MiscUtil.PrintError("Failed to setup a writer for output fie %s " % OptionsInfo["OutfileFlexRes"])
801
802 return (Writer, WriterFlexRes)
803
804
805 def CloseWriters(Writer, WriterFlexRes):
806 """Close writers."""
807
808 if Writer is not None:
809 Writer.close()
810
811 if WriterFlexRes is not None:
812 WriterFlexRes.close()
813
814
815 def ComputeGridCenter(LigandFile):
816 """Compute grid center from ligand file."""
817
818 GridCenter = []
819
820 MiscUtil.PrintInfo("\nComputing grid center from reference ligand file %s..." % LigandFile)
821
822 Mols = RDKitUtil.ReadMolecules(LigandFile)
823 Mol = Mols[0]
824
825 Centroid = rdMolTransforms.ComputeCentroid(Mol.GetConformer())
826
827 GridCenter = [Centroid.x, Centroid.y, Centroid.z]
828
829 GridCenterFormatted = ["%.3f" % Value for Value in GridCenter]
830 MiscUtil.PrintInfo("GridCenter: %s" % (" ".join(GridCenterFormatted)))
831
832 return GridCenter
833
834
835 def ProcessReceptorOptions():
836 """Process receptor options."""
837
838 ReceptorFile, ReceptorMapsPrefix, UseReceptorFile, UseReceptorMapsPrefix = [None, None, False, False]
839
840 Receptor = Options["--receptor"]
841 if os.path.isfile(Receptor):
842 ReceptorFile = Receptor
843 UseReceptorFile = True
844 else:
845 ReceptorMapsPrefix = Receptor
846 UseReceptorMapsPrefix = True
847
848 OptionsInfo["Receptor"] = Receptor
849
850 OptionsInfo["ReceptorFile"] = ReceptorFile
851 OptionsInfo["UseReceptorFile"] = UseReceptorFile
852
853 OptionsInfo["ReceptorMapsPrefix"] = ReceptorMapsPrefix
854 OptionsInfo["UseReceptorMapsPrefix"] = UseReceptorMapsPrefix
855
856 UseReceptorFlexFile = False
857 ReceptorFlexFile = None
858 if not re.match("^None$", Options["--receptorFlexFile"], re.I):
859 UseReceptorFlexFile = True
860 ReceptorFlexFile = Options["--receptorFlexFile"]
861 OptionsInfo["ReceptorFlexFile"] = ReceptorFlexFile
862 OptionsInfo["UseReceptorFlexFile"] = UseReceptorFlexFile
863
864 if OptionsInfo["UseAD4Forcefield"]:
865 if OptionsInfo["UseReceptorFile"]:
866 MiscUtil.PrintError(
867 'The value specified, %s, for option "-r, --receptor" is not valid for %s forcefield. Supported value: Affinity maps prefix.'
868 % (OptionsInfo["Receptor"], Options["--forcefield"])
869 )
870
871
872 def ProcessForcefieldWeightParamatersOption(ParamsOptionName, ParamsOptionValue, ParamsDefaultInfo=None):
873 """Process forcefield weight paramaters option."""
874
875 ParamsInfo = {
876 "AD4Vdw": 0.1662,
877 "AD4HydrogenBond": 0.1209,
878 "AD4Electrostatic": 0.1406,
879 "AD4Desolvation": 0.1322,
880 "AD4GlueLinearAttraction": 50.0,
881 "AD4Rot": 0.2983,
882 "VinaGaussian1": -0.035579,
883 "VinaGaussian2": -0.005156,
884 "VinaRepulsion": 0.840245,
885 "VinaHydrophobic": -0.035069,
886 "VinaHydrogenBond": -0.587439,
887 "VinaGlueLinearAttraction": 50.0,
888 "VinaRot": 0.05846,
889 "VinardoGaussian1": -0.045,
890 "VinardoRepulsion": 0.8,
891 "VinardoHydrophobic": -0.035,
892 "VinardoHydrogenBond": -0.600,
893 "VinardoGlueLinearAttraction": 50.0,
894 "VinardoRot": 0.05846,
895 }
896
897 # Setup a canonical paramater names...
898 ValidParamNames = []
899 CanonicalParamNamesMap = {}
900 for ParamName in sorted(ParamsInfo):
901 ValidParamNames.append(ParamName)
902 CanonicalParamNamesMap[ParamName.lower()] = ParamName
903
904 # Update default values...
905 if ParamsDefaultInfo is not None:
906 for ParamName in ParamsDefaultInfo:
907 if ParamName not in ParamsInfo:
908 MiscUtil.PrintError(
909 'The default parameter name, %s, specified using "%s" to function ProcessForcefieldWeightParamatersOption is not a valid name. Supported parameter names: %s'
910 % (ParamName, ParamsDefaultInfo, " ".join(ValidParamNames))
911 )
912 ParamsInfo[ParamName] = ParamsDefaultInfo[ParamName]
913
914 if re.match("^auto$", ParamsOptionValue, re.I):
915 return ParamsInfo
916
917 ParamsOptionValueWords = ParamsOptionValue.split(",")
918 if len(ParamsOptionValueWords) % 2:
919 MiscUtil.PrintError(
920 'The number of comma delimited paramater names and values, %d, specified using "%s" option must be an even number.'
921 % (len(ParamsOptionValueWords), ParamsOptionName)
922 )
923
924 # Validate paramater name and value pairs...
925 for Index in range(0, len(ParamsOptionValueWords), 2):
926 Name = ParamsOptionValueWords[Index].strip()
927 Value = ParamsOptionValueWords[Index + 1].strip()
928
929 CanonicalName = Name.lower()
930 if CanonicalName not in CanonicalParamNamesMap:
931 MiscUtil.PrintError(
932 'The parameter name, %s, specified using "%s" is not a valid name. Supported parameter names: %s'
933 % (Name, ParamsOptionName, " ".join(ValidParamNames))
934 )
935
936 ParamName = CanonicalParamNamesMap[CanonicalName]
937
938 if not MiscUtil.IsFloat(Value):
939 MiscUtil.PrintError(
940 'The parameter value, %s, specified for parameter name, %s, using "%s" must be a float.'
941 % (Value, Name, ParamsOptionName)
942 )
943 ParamValue = float(Value)
944
945 ParamsInfo[ParamName] = ParamValue
946
947 return ParamsInfo
948
949
950 def ProcessModeOption():
951 """Process mode option."""
952
953 Mode = Options["--mode"]
954 DockMode, LocalOptimizationOnlyMode, ScoreOnlyMode = [False] * 3
955 if re.match("^Dock$", Mode, re.I):
956 Mode = "Dock"
957 DockMode = True
958 elif re.match("^LocalOptimizationOnly$", Mode, re.I):
959 Mode = "LocalOptimizationOnly"
960 LocalOptimizationOnlyMode = True
961 elif re.match("^ScoreOnly$", Mode, re.I):
962 Mode = "ScoreOnly"
963 ScoreOnlyMode = True
964 else:
965 MiscUtil.PrintError(
966 'The value specified, %s, for option "-m, --mode" is not valid. Supported values: Dock LocalOptimizationOnly ScoreOnly'
967 % OptionsInfo["Mode"]
968 )
969
970 OptionsInfo["Mode"] = Mode
971 OptionsInfo["DockMode"] = DockMode
972 OptionsInfo["LocalOptimizationOnlyMode"] = LocalOptimizationOnlyMode
973 OptionsInfo["ScoreOnlyMode"] = ScoreOnlyMode
974
975
976 def ProcessEnergyLabelOptions():
977 """Process energy label options."""
978
979 Forcefield = OptionsInfo["Forcefield"]
980
981 EnergyLabel = Options["--energyLabel"]
982 if re.match("^auto$", EnergyLabel, re.I):
983 EnergyLabel = "%s_Total_Energy (kcal/mol)" % Forcefield
984 OptionsInfo["EnergyLabel"] = EnergyLabel
985
986 EnergyLabelsList = []
987 DefaultLabels = [
988 "%s_Intermolecular_Energy (kcal/mol)" % (Forcefield),
989 "%s_Internal_Energy (kcal/mol)" % (Forcefield),
990 "%s_Torsions_Energy (kcal/mol)" % (Forcefield),
991 ]
992
993 EnergyLabels = Options["--energyComponentsLabels"]
994 if not re.match("^auto$", EnergyLabels, re.I):
995 EnergyLabelsWords = EnergyLabels.split(",")
996 if len(EnergyLabelsWords) != 3:
997 MiscUtil.PrintError(
998 'The specified value, %s, for option "--energyComponentsLabels " is not valid. It must contain 3 text values separated by comma.'
999 % EnergyLabels
1000 )
1001
1002 for LabelIndex, Label in enumerate(EnergyLabelsWords):
1003 Label = Label.strip()
1004 if re.match("^auto$", Label, re.I):
1005 Label = DefaultLabels[LabelIndex]
1006
1007 EnergyLabelsList.append(Label)
1008 else:
1009 EnergyLabelsList = DefaultLabels
1010
1011 OptionsInfo["EnergyComponentsLabels"] = EnergyLabels
1012 OptionsInfo["EnergyComponentsLabelsList"] = EnergyLabelsList
1013
1014 SetupEnergyLabelsAndIndices()
1015
1016
1017 def SetupEnergyLabelsAndIndices():
1018 """Setup energy data labels and indices for writing out energy values."""
1019
1020 EnergyLabelsList = []
1021 EnergyValueIndicesList = []
1022
1023 # Total energy is always at index 0 in the list retured by vina.score (),
1024 # vina.optimize(), and vina.energies() during ScoreOnly, LocalOptimizationOnly
1025 # and Dock modes...
1026 EnergyLabelsList.append(OptionsInfo["EnergyLabel"])
1027 EnergyValueIndicesList.append(0)
1028
1029 OptionsInfo["EnergyLabelsList"] = EnergyLabelsList
1030 OptionsInfo["EnergyValueIndicesList"] = EnergyValueIndicesList
1031
1032 if not OptionsInfo["EnergyComponents"]:
1033 return
1034
1035 # Setup energy labels and indices for energy components...
1036 if OptionsInfo["ScoreOnlyMode"]:
1037 # vina.score returns a list containing the following values:
1038 #
1039 # Vina/Vinardo FF: columns=[total, lig_inter, flex_inter, other_inter, flex_intra, lig_intra, torsions, lig_intra best pose]
1040 # AutoDock FF: [total, lig_inter, flex_inter, other_inter, flex_intra, lig_intra, torsions, -lig_intra]
1041 #
1042 IntermolecularEnergyIndex, IntramolecularEnergyIndex, TorsionEnergyIndex = [1, 5, 6]
1043 elif OptionsInfo["LocalOptimizationOnlyMode"]:
1044 # vina.optimize returns a list containing the following values:
1045 #
1046 # Vina/Vinardo FF: [total, lig_inter, flex_inter, other_inter, flex_intra, lig_intra, torsions, lig_intra best pose]
1047 # AutoDock FF: [total, lig_inter, flex_inter, other_inter, flex_intra, lig_intra, torsions, -lig_intra]
1048 #
1049 IntermolecularEnergyIndex, IntramolecularEnergyIndex, TorsionEnergyIndex = [1, 5, 6]
1050 elif OptionsInfo["DockMode"]:
1051 # vina.energies returns a list containing the following values:
1052 #
1053 # Vina/Vinardo FF: [total, inter, intra, torsions, intra best pose]
1054 # AutoDock FF: [total, inter, intra, torsions, -intra]
1055 #
1056 IntermolecularEnergyIndex, IntramolecularEnergyIndex, TorsionEnergyIndex = [1, 2, 3]
1057 else:
1058 MiscUtil.PrintError(
1059 'The value specified, %s, for option "-m, --mode" is not valid. Supported values: Dock LocalOptimizationOnly ScoreOnly'
1060 % OptionsInfo["Mode"]
1061 )
1062
1063 OptionsInfo["EnergyLabelsList"].extend(OptionsInfo["EnergyComponentsLabelsList"])
1064 OptionsInfo["EnergyValueIndicesList"].extend(
1065 [IntermolecularEnergyIndex, IntramolecularEnergyIndex, TorsionEnergyIndex]
1066 )
1067
1068
1069 def ProcessOptions():
1070 """Process and validate command line arguments and options."""
1071
1072 MiscUtil.PrintInfo("Processing options...")
1073
1074 # Validate options...
1075 ValidateOptions()
1076
1077 OptionsInfo["Infile"] = Options["--infile"]
1078 ParamsDefaultInfoOverride = {"RemoveHydrogens": False}
1079 OptionsInfo["InfileParams"] = MiscUtil.ProcessOptionInfileParameters(
1080 "--infileParams",
1081 Options["--infileParams"],
1082 InfileName=Options["--infile"],
1083 ParamsDefaultInfo=ParamsDefaultInfoOverride,
1084 )
1085
1086 OptionsInfo["Outfile"] = Options["--outfile"]
1087 OptionsInfo["OutfileParams"] = MiscUtil.ProcessOptionOutfileParameters(
1088 "--outfileParams", Options["--outfileParams"]
1089 )
1090
1091 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--outfile"])
1092 OptionsInfo["OutfileFlexRes"] = "%s_Flex_Receptor.%s" % (FileName, FileExt)
1093
1094 GridCenterLigandFile, GridCenterList, GridCenterByLigandFile = [None, None, False]
1095 GridCenter = Options["--gridCenter"]
1096 if re.search(",", GridCenter, re.I):
1097 GridCenterList = [float(Value) for Value in GridCenter.split(",")]
1098 else:
1099 GridCenterByLigandFile = True
1100 GridCenterLigandFile = GridCenter
1101 GridCenterList = ComputeGridCenter(GridCenterLigandFile)
1102
1103 OptionsInfo["GridCenter"] = GridCenter
1104 OptionsInfo["GridCenterList"] = GridCenterList
1105 OptionsInfo["GridCenterLigandFile"] = GridCenterLigandFile
1106 OptionsInfo["GridCenterByLigandFile"] = GridCenterByLigandFile
1107
1108 OptionsInfo["EnergyComponents"] = True if re.match("^yes$", Options["--energyComponents"], re.I) else False
1109
1110 OptionsInfo["EnergyRange"] = float(Options["--energyRange"])
1111 OptionsInfo["Exhaustiveness"] = int(Options["--exhaustiveness"])
1112
1113 Forcefield = Options["--forcefield"]
1114 UseAD4Forcefield, UseVinaForcefield, UseVinardoForcefield = [False, False, False]
1115 if re.match("^AD4$", Forcefield, re.I):
1116 Forcefield = "AD4"
1117 UseAD4Forcefield = True
1118 elif re.match("^Vina$", Forcefield, re.I):
1119 Forcefield = "Vina"
1120 UseVinaForcefield = True
1121 elif re.match("^Vinardo$", Forcefield, re.I):
1122 Forcefield = "Vinardo"
1123 UseVinardoForcefield = True
1124 else:
1125 MiscUtil.PrintError(
1126 'The value specified, %s, for option "-f, --forcefield" is not valid. Supported values: AD4 Vina Vinardo'
1127 )
1128 OptionsInfo["Forcefield"] = Forcefield
1129 OptionsInfo["UseAD4Forcefield"] = UseAD4Forcefield
1130 OptionsInfo["UseVinaForcefield"] = UseVinaForcefield
1131 OptionsInfo["UseVinardoForcefield"] = UseVinardoForcefield
1132
1133 OptionsInfo["ForcefieldWeightParams"] = ProcessForcefieldWeightParamatersOption(
1134 "--forcefieldWeightParams", Options["--forcefieldWeightParams"]
1135 )
1136 OptionsInfo["ForcefieldWeightParamsSpecified"] = (
1137 False if re.match("^auto$", Options["--forcefieldWeightParams"], re.I) else True
1138 )
1139
1140 GridSize = Options["--gridSize"]
1141 GridSizeList = [float(Value) for Value in GridSize.split(",")]
1142 OptionsInfo["GridSize"] = GridSize
1143 OptionsInfo["GridSizeList"] = GridSizeList
1144
1145 OptionsInfo["GridSpacing"] = float(Options["--gridSpacing"])
1146
1147 OptionsInfo["MaxEvaluations"] = int(Options["--maxEvaluations"])
1148 OptionsInfo["MinRMSD"] = float(Options["--minRMSD"])
1149
1150 MergeHydrogens = Options["--mergeHydrogens"]
1151 if re.match("^auto$", MergeHydrogens, re.I):
1152 MergeHydrogens = True if OptionsInfo["UseAD4Forcefield"] else False
1153 else:
1154 MergeHydrogens = True if re.match("^yes$", MergeHydrogens, re.I) else False
1155 OptionsInfo["MergeHydrogens"] = MergeHydrogens
1156
1157 ProcessModeOption()
1158
1159 OptionsInfo["MPMode"] = True if re.match("^yes$", Options["--mp"], re.I) else False
1160 OptionsInfo["MPParams"] = MiscUtil.ProcessOptionMultiprocessingParameters("--mpParams", Options["--mpParams"])
1161
1162 OptionsInfo["NumPoses"] = int(Options["--numPoses"])
1163
1164 if re.match("^auto$", Options["--numThreads"], re.I):
1165 NumThreads = 1 if OptionsInfo["MPMode"] else 0
1166 else:
1167 NumThreads = int(Options["--numThreads"])
1168 OptionsInfo["NumThreads"] = NumThreads
1169
1170 OptionsInfo["Precision"] = int(Options["--precision"])
1171 OptionsInfo["RandomSeed"] = int(Options["--randomSeed"])
1172
1173 OptionsInfo["SkipRefinement"] = True if re.match("^yes$", Options["--skipRefinement"], re.I) else False
1174
1175 OptionsInfo["ValidateMolecules"] = True if re.match("^yes$", Options["--validateMolecules"], re.I) else False
1176
1177 if re.match("^auto$", Options["--vinaVerbosity"], re.I):
1178 VinaVerbosity = 0 if OptionsInfo["MPMode"] else 1
1179 else:
1180 VinaVerbosity = int(Options["--vinaVerbosity"])
1181 OptionsInfo["VinaVerbosity"] = VinaVerbosity
1182
1183 OptionsInfo["Overwrite"] = Options["--overwrite"]
1184 OptionsInfo["QuietMode"] = True if re.match("^yes$", Options["--quiet"], re.I) else False
1185
1186 ProcessEnergyLabelOptions()
1187 ProcessReceptorOptions()
1188
1189
1190 def RetrieveOptions():
1191 """Retrieve command line arguments and options."""
1192
1193 # Get options...
1194 global Options
1195 Options = docopt(_docoptUsage_)
1196
1197 # Set current working directory to the specified directory...
1198 WorkingDir = Options["--workingdir"]
1199 if WorkingDir:
1200 os.chdir(WorkingDir)
1201
1202 # Handle examples option...
1203 if "--examples" in Options and Options["--examples"]:
1204 MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
1205 sys.exit(0)
1206
1207
1208 def ValidateOptions():
1209 """Validate option values."""
1210
1211 MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
1212 MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "sdf sd mol")
1213
1214 MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "sdf sd")
1215 MiscUtil.ValidateOptionsOutputFileOverwrite(
1216 "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
1217 )
1218 MiscUtil.ValidateOptionsDistinctFileNames(
1219 "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
1220 )
1221
1222 FileDir, FileName, FileExt = MiscUtil.ParseFileName(Options["--receptor"])
1223 if not MiscUtil.IsEmpty(FileExt):
1224 MiscUtil.ValidateOptionFilePath("-r, --receptor", Options["--receptor"])
1225 MiscUtil.ValidateOptionFileExt("-r, --receptor", Options["--receptor"], "pdbqt")
1226 else:
1227 AffinityMapFiles = glob.glob("%s.*.map" % Options["--receptor"])
1228 if len(AffinityMapFiles) == 0:
1229 MiscUtil.PrintError(
1230 'The receptor affinity map files, %s.*.map, corresponding to maps prefix, %s, specified using option, "-r, --receptor" option don\'t exist.'
1231 % (Options["--receptor"], Options["--receptor"])
1232 )
1233
1234 if not re.match("^None$", Options["--receptorFlexFile"], re.I):
1235 MiscUtil.ValidateOptionFilePath("--receptorFlexFile", Options["--receptorFlexFile"])
1236 MiscUtil.ValidateOptionFileExt("--receptorFlexFile", Options["--receptorFlexFile"], "pdbqt")
1237 if os.path.isfile(Options["--receptor"]):
1238 MiscUtil.ValidateOptionsDistinctFileNames(
1239 "-r, --receptor", Options["--receptor"], "--receptorFlexFile", Options["--receptorFlexFile"]
1240 )
1241
1242 if re.search(",", Options["--gridCenter"], re.I):
1243 MiscUtil.ValidateOptionNumberValues("-g, --gridCenter", Options["--gridCenter"], 3, ",", "float", {">=": 0.0})
1244 else:
1245 MiscUtil.ValidateOptionFilePath("-g, --gridCenter", Options["--gridCenter"])
1246 MiscUtil.ValidateOptionFileExt("-g, --gridCenter", Options["--gridCenter"], "sdf sd mol pdb")
1247
1248 MiscUtil.ValidateOptionTextValue("--energyComponents", Options["--energyComponents"], "yes no")
1249
1250 MiscUtil.ValidateOptionFloatValue("--energyRange", Options["--energyRange"], {">": 0})
1251
1252 MiscUtil.ValidateOptionIntegerValue("--exhaustiveness", Options["--exhaustiveness"], {">": 0})
1253
1254 MiscUtil.ValidateOptionTextValue("-f, --forcefield", Options["--forcefield"], "AD4 Vina Vinardo")
1255
1256 MiscUtil.ValidateOptionNumberValues("--gridSize", Options["--gridSize"], 3, ",", "float", {">": 0.0})
1257 MiscUtil.ValidateOptionFloatValue("--gridSpacing", Options["--gridSpacing"], {">": 0})
1258
1259 MiscUtil.ValidateOptionIntegerValue("--maxEvaluations", Options["--maxEvaluations"], {">=": 0})
1260 MiscUtil.ValidateOptionFloatValue("--minRMSD", Options["--minRMSD"], {">": 0})
1261
1262 MiscUtil.ValidateOptionTextValue("--mergeHydrogens", Options["--mergeHydrogens"], "yes no auto")
1263
1264 MiscUtil.ValidateOptionTextValue("-m, --mode", Options["--mode"], "Dock LocalOptimizationOnly ScoreOnly")
1265
1266 MiscUtil.ValidateOptionTextValue("--mp", Options["--mp"], "yes no")
1267
1268 MiscUtil.ValidateOptionIntegerValue("--numPoses", Options["--numPoses"], {">": 0})
1269
1270 if not re.match("^auto$", Options["--numThreads"], re.I):
1271 MiscUtil.ValidateOptionIntegerValue("--numThreads", Options["--numThreads"], {">": 0})
1272
1273 MiscUtil.ValidateOptionIntegerValue("-p, --precision", Options["--precision"], {">": 0})
1274 MiscUtil.ValidateOptionIntegerValue("--randomSeed", Options["--randomSeed"], {})
1275
1276 MiscUtil.ValidateOptionTextValue("--skipRefinement", Options["--skipRefinement"], "yes no")
1277
1278 MiscUtil.ValidateOptionTextValue("-v, --validateMolecules", Options["--validateMolecules"], "yes no")
1279
1280 if not re.match("^auto$", Options["--vinaVerbosity"], re.I):
1281 MiscUtil.ValidateOptionIntegerValue("--vinaVerbosity", Options["--vinaVerbosity"], {">=": 0})
1282 MiscUtil.ValidateOptionTextValue("--vinaVerbosity", Options["--vinaVerbosity"], "0 1 2")
1283
1284 MiscUtil.ValidateOptionTextValue("-q, --quiet", Options["--quiet"], "yes no")
1285
1286
1287 # Setup a usage string for docopt...
1288 _docoptUsage_ = """
1289 VinaPerformDocking.py - Perform docking
1290
1291 Usage:
1292 VinaPerformDocking.py [--energyComponents <yes or no>] [--energyComponentsLabels <Label1, label2, Label3>]
1293 [--energyLabel <text>] [--energyRange <number>] [--exhaustiveness <number>]
1294 [--forcefield <AD4, Vina, or Vinardo>] [--forcefieldWeightParams <Name,Value,...>] [--gridSpacing <number>]
1295 [--gridSize <xsize,ysize,zsize>] [--infileParams <Name,Value,...>] [--maxEvaluations <number>]
1296 [--mergeHydrogens <yes or no>] [--minRMSD <number>] [--mode <Dock, LocalOptimizationOnly, ScoreOnly>] [--mp <yes or no>]
1297 [--mpParams <Name, Value,...>] [--numPoses <number>] [--numThreads <number>] [--outfileParams <Name,Value,...> ]
1298 [--overwrite] [--precision <number>] [--quiet <yes or no>] [--randomSeed <number>] [--receptorFlexFile <receptor flex file>]
1299 [--skipRefinement <yes or no>] [--validateMolecules <yes or no>] [--vinaVerbosity <number>]
1300 [-w <dir>] -g <RefLigandFile or x,y,z> -r <receptorfile or maps prerfix> -i <infile> -o <outfile>
1301 VinaPerformDocking.py -h | --help | -e | --examples
1302
1303 Description:
1304 Dock molecules against a protein receptor using AutoDock Vina [Ref 168, 169].
1305 The molecules must have 3D coordinates in input file. In addition, the hydrogens
1306 must be present for all molecules in input file.
1307
1308 No protein receptor preparation is performed during docking. It must be prepared
1309 employing standalone scripts available as part of AutoDock Vina. You may
1310 optionally specify flexible residues in the binding pocket to prepare a flexible
1311 receptor file and employ it for docking molecules along with the fixed receptor
1312 file.
1313
1314 The following three forecefileds are available to score molecules: AD4
1315 (AutoDock4), Vina, and Vinardo (Vina RaDii Optimized) [Ref 170].
1316
1317 The supported input file formats are shown below:
1318
1319 Rigid/Flexible protein receptor files - PDBQT(.pdbqt)
1320 Reference ligand file: PDB(.pdb), Mol (.mol), SD (.sdf, .sd)
1321
1322 Input molecules file - Mol (.mol), SD (.sdf, .sd)
1323
1324 The supported output file format is: SD (.sdf, .sd).
1325
1326 The following output files are generated:
1327
1328 <OutfileRoot>.<OutfileExt> - Docked/scored molecules
1329 <OutfileRoot>_Flex_Receptor.<OutfileExt> - Docked poses for flexible
1330 residues
1331
1332 The flexible receptor output file contains docked poses corresponding to
1333 flexible residues. It is only generated during 'Dock' value of '-m, --mode'
1334 option. The number of poses in this file matches those written to the
1335 output file containing docked molecules.
1336
1337 Options:
1338 --energyComponents <yes or no> [default: no]
1339 Write out binding energy components of the total binding energy docking
1340 score to outfile. The following three energy components are written to
1341 outfile: intermolecula energy, internal energy, and torsions energy.
1342 --energyComponentsLabels <Label1, label2, Label3> [default: auto]
1343 A triplet of comma delimited values corresponding to energy data field
1344 labels for writing out the binding energy components to outfile. You must
1345 specify all three values. A value of 'None' implies the use of the default
1346 labels as shown below:
1347
1348 Label1: <ForcefieldName>_Intermolecular_Energy (kcal/mol)
1349 Label2: <ForcefieldName>_Internal_Energy (kcal/mol)
1350 Label3: <ForcefieldName>_Torsions_Energy (kcal/mol)
1351
1352 --energyLabel <text> [default: auto]
1353 Energy data field label for writing out binding energy docking score to
1354 output file. Default: <ForcefieldName>_Total_Energy (kcal/mol).
1355 --energyRange <number> [default: 3.0]
1356 Maximum energy difference from the best pose during the generation of
1357 poses. Units: kcal/mol.
1358 -e, --examples
1359 Print examples.
1360 --exhaustiveness <number> [default: 8]
1361 Exhaustiveness of global MC search. The higher values make the search
1362 more exhaustive and it takes longer to complete. You may want to use
1363 '16' or '32' as the value of '--exhaustiveness ' to increase the accuracy of
1364 your pose prediction.
1365 -f, --forcefield <AD4, Vina, or Vinardo> [default: Vina]
1366 Forcefield to use for scoring. Possible values: AD4 (AutoDock 4), Vina
1367 [Ref 169, 169], or Vinardo (Vina RaDii Optimized) [Ref 170].
1368
1369 You must specify affinity maps using '-r, --receptor' option during the use
1370 of 'AD4' forcefield.
1371 --forcefieldWeightParams <Name,Value,...> [default: auto]
1372 A comma delimited list of parameter name and value pairs for forcefield
1373 scoring.
1374
1375 The supported parameter names along with their default values are
1376 are shown below for different forcefields:
1377
1378 AD4 (6 weights):
1379
1380 ad4Vdw, 0.1662, ad4HydrogenBond, 0.1209, ad4Electrostatic, 0.1406,
1381 ad4Desolvation, 0.1322, ad4GlueLinearAttraction, 50.0,
1382 ad4Rot, 0.2983
1383
1384 Vina (7 weights):
1385
1386 vinaGaussian1, -0.035579, vinaGaussian2, -0.005156,
1387 vinaRepulsion, 0.840245, vinaHydrophobic, -0.035069,
1388 vinaHydrogenBond, -0.587439, vinaGlueLinearAttraction, 50.0,
1389 vinaRot, 0.05846
1390
1391 Vinardo (6 weights):
1392
1393 vinardoGaussian1, -0.045, vinardoRepulsion, 0.8,
1394 vinardoHydrophobic, -0.035, vinardoHydrogenBond, -0.600
1395 vinardoGlueLinearAttraction, 50.0, vinardoRot, 0.05846
1396
1397 The glue weight parameter corresponds to linear attraction for macrocycle
1398 closure and has the same value for AD4, Vina, and Vinardo. The rot weight
1399 has the same value for Vina and Vinardo.
1400 -g, --gridCenter <RefLigandFile or x,y,z>
1401 Reference ligand file for calculating the docking grid center or a triplet
1402 of comma delimited values in Angstrom corresponding to grid center.
1403
1404 This is required option. However, it is ignored during the specification
1405 of maps prefix for '-r, --receptor' option.
1406 --gridSize <xsize,ysize,zsize> [default: 25.0, 25.0, 25.0]
1407 Docking grid size in Angstrom.
1408 --gridSpacing <number> [default: 0.375]
1409 Docking grid spacing in Angstrom.
1410 -h, --help
1411 Print this help message.
1412 -i, --infile <infile>
1413 Input file name containing molecules for docking against a receptor. The
1414 molecules must have 3D coordinates in input file. In addition, the hydrogens
1415 must be present for all molecules in input file. The input file may contain 3D
1416 conformers.
1417 --infileParams <Name,Value,...> [default: auto]
1418 A comma delimited list of parameter name and value pairs for reading
1419 molecules from files. The supported parameter names for different file
1420 formats, along with their default values, are shown below:
1421
1422 SD, MOL: removeHydrogens,no,sanitize,yes,strictParsing,yes
1423
1424 --maxEvaluations <number> [default: 0]
1425 Maximum number of evaluations to perform for each MC run during docking.
1426 By default, its value is set 0 and the number of MC evaluations is determined
1427 using heuristic rules.
1428 --mergeHydrogens <yes or no> [default: auto]
1429 Merge hydrogens during preparation of molecules for docking. The hydrogens
1430 are automatically merged during 'AD4' value of '-f, --forcefield' option and its
1431 value is set to 'yes'. Otherwise, it's set to 'no'.
1432 --minRMSD <number> [default: 1.0]
1433 Minimum RMSD between output poses in Angstrom.
1434 -m, --mode <Dock, LocalOptimizationOnly, ScoreOnly> [default: Dock]
1435 Dock molecules or simply score molecules without performing any docking.
1436 The supported values along with a brief explanation of the expected
1437 behavior are shown below:
1438
1439 Dock: Global search along with local optimization and scoring after
1440 docking
1441 LocalOptimizationOnly: Local optimization and scoring without any docking
1442 ScoreOnly: Scoring without any local optimizatiom and docking
1443
1444 The 'ScoreOnly" allows you to score 3D moleculed from input file which
1445 are already positioned in a binding pocket of a receptor.
1446 --mp <yes or no> [default: no]
1447 Use multiprocessing.
1448
1449 By default, input data is retrieved in a lazy manner via mp.Pool.imap()
1450 function employing lazy RDKit data iterable. This allows processing of
1451 arbitrary large data sets without any additional requirements memory.
1452
1453 All input data may be optionally loaded into memory by mp.Pool.map()
1454 before starting worker processes in a process pool by setting the value
1455 of 'inputDataMode' to 'InMemory' in '--mpParams' option.
1456
1457 A word to the wise: The default 'chunkSize' value of 1 during 'Lazy' input
1458 data mode may adversely impact the performance. The '--mpParams' section
1459 provides additional information to tune the value of 'chunkSize'.
1460 --mpParams <Name,Value,...> [default: auto]
1461 A comma delimited list of parameter name and value pairs to configure
1462 multiprocessing.
1463
1464 The supported parameter names along with their default and possible
1465 values are shown below:
1466
1467 chunkSize, auto
1468 inputDataMode, Lazy [ Possible values: InMemory or Lazy ]
1469 numProcesses, auto [ Default: mp.cpu_count() ]
1470
1471 These parameters are used by the following functions to configure and
1472 control the behavior of multiprocessing: mp.Pool(), mp.Pool.map(), and
1473 mp.Pool.imap().
1474
1475 The chunkSize determines chunks of input data passed to each worker
1476 process in a process pool by mp.Pool.map() and mp.Pool.imap() functions.
1477 The default value of chunkSize is dependent on the value of 'inputDataMode'.
1478
1479 The mp.Pool.map() function, invoked during 'InMemory' input data mode,
1480 automatically converts RDKit data iterable into a list, loads all data into
1481 memory, and calculates the default chunkSize using the following method
1482 as shown in its code:
1483
1484 chunkSize, extra = divmod(len(dataIterable), len(numProcesses) * 4)
1485 if extra: chunkSize += 1
1486
1487 For example, the default chunkSize will be 7 for a pool of 4 worker processes
1488 and 100 data items.
1489
1490 The mp.Pool.imap() function, invoked during 'Lazy' input data mode, employs
1491 'lazy' RDKit data iterable to retrieve data as needed, without loading all the
1492 data into memory. Consequently, the size of input data is not known a priori.
1493 It's not possible to estimate an optimal value for the chunkSize. The default
1494 chunkSize is set to 1.
1495
1496 The default value for the chunkSize during 'Lazy' data mode may adversely
1497 impact the performance due to the overhead associated with exchanging
1498 small chunks of data. It is generally a good idea to explicitly set chunkSize to
1499 a larger value during 'Lazy' input data mode, based on the size of your input
1500 data and number of processes in the process pool.
1501
1502 The mp.Pool.map() function waits for all worker processes to process all
1503 the data and return the results. The mp.Pool.imap() function, however,
1504 returns the the results obtained from worker processes as soon as the
1505 results become available for specified chunks of data.
1506
1507 The order of data in the results returned by both mp.Pool.map() and
1508 mp.Pool.imap() functions always corresponds to the input data.
1509 -n, --numPoses <number> [default: 1]
1510 Number of docked poses to generate for each molecule and write out
1511 to output file. This option is only valid for 'Dock' value of '-m, --mode'
1512 option.
1513 --numThreads <number> [default: auto]
1514 Number of threads/CPUs to use for MC search calculation in Vina. The
1515 default value is set to 1 during multiprocessing for 'Yes' value of '--mp'
1516 option. Otherwise, it's set to 0 and rely on Vina detect and use all
1517 available CPUs for multi-threading.
1518 -o, --outfile <outfile>
1519 Output file name for writing out molecules. The flexible receptor residues
1520 are written to <OutfileRoot>_Flex_Receptor.<OutfileExt>.
1521 --outfileParams <Name,Value,...> [default: auto]
1522 A comma delimited list of parameter name and value pairs for writing
1523 molecules to files. The supported parameter names for different file
1524 formats, along with their default values, are shown below:
1525
1526 SD: kekulize,yes,forceV3000,no
1527
1528 --overwrite
1529 Overwrite existing files.
1530 --precision <number> [default: 2]
1531 Floating point precision for writing energy values.
1532 -q, --quiet <yes or no> [default: no]
1533 Use quiet mode. The warning and information messages will not be printed.
1534 --randomSeed <number> [default: 0]
1535 Random seed for MC search calculations. A value of zero implies it's randomly
1536 chosen by Vina during the calculation.
1537 -r, --receptor <receptor file or maps prerfix>
1538 Protein receptor file name or prefix for affinity map files corresponding
1539 to the fixed portion of the receptor.
1540
1541 You must specify affinity map files for 'AD4' forcefield. The affinity map
1542 files correspond to <MapsPrefix>.*.map and must be present
1543
1544 The supported receptor file format is PDBQT (.pdbqt). It must contain a
1545 prepared protein receptor ready for docking. You may prepare a PDBQT
1546 receptor file from a PDB file employing the command line scripts available
1547 with AutoDock Vina and Meeko. For example: prepare_receptor,
1548 prepare_flexreceptor.py, or or mk_make_recepror.py.
1549
1550 You may want to perform the following steps to clean up your PDB file
1551 before generating a PDBQT receptor file: Remove extraneous molecules
1552 such as solvents, ions, and ligand etc.; Extract data for a chain containing
1553 the binding pocket of interest.
1554 --receptorFlexFile <receptor flex file> [default: none]
1555 Protein receptor file name corresponding to the flexible portion of the
1556 receptor. The supported receptor file format is PDBQT (.pdbqt). It must
1557 contain a prepared protein receptor ready for docking. You may prepare
1558 a flexible PDBQT receptor file from a PDB file employing the command line
1559 script prepare_flexreceptor or mk_make_recepror.py available with Autodock
1560 Vina and Meeko.
1561 --skipRefinement <yes or no> [default: no]
1562 Skip refinement. Vina is initialized to skip the use of explicit receptor atoms,
1563 instead of precalculated grids, during the following three docking modes:
1564 Dock, LocalOptimizationOnly, and ScoreOnly.
1565 -v, --validateMolecules <yes or no> [default: yes]
1566 Validate molecules for docking. The input molecules must have 3D coordinates
1567 and the hydrogens must be present. You may skip validation of molecules in
1568 input file containing all valid molecules.
1569 --vinaVerbosity <number> [default: auto]
1570 Verbosity level for running Vina. Possible values: 0 - No output; 1 - Normal
1571 output; 2 - Verbose. Default: 0 during multiprocessing for 'Yes' value of '--mp'
1572 option; otherwise, its set to 1. A non-zero value is not recommended during
1573 multiprocessing. It doesn't work due to the mingling of the Vina output
1574 from multiple processes.
1575 -w, --workingdir <dir>
1576 Location of working directory which defaults to the current directory.
1577
1578 Examples:
1579 To dock molecules using Vina forcefield against a prepared receptor
1580 corresponding to chain A extracted from PDB file 1R4L.pdb, multi-threading
1581 across all available CPUs during Vina MC search, calculating grid center form
1582 a reference ligand file, using grid box size of 25 Angstrom, generating one
1583 pose for each molecule, and write out a SD file containing docked molecules,
1584 type:
1585
1586 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1587 -r SampleACE2Receptor.pdbqt -i SampleACE2Ligands.sdf
1588 -o SampleACE2LigandsOut.sdf
1589
1590 To run the first example for generating multiple docked poses for each
1591 molecule and write out all energy terms to a SD file, type:
1592
1593 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1594 -r SampleACE2Receptor.pdbqt --numPoses 5 --energyComponents yes
1595 -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1596
1597 To run the first example for docking molecules using Vinardo forcefield and write
1598 out a SD file, type:
1599
1600 % VinaPerformDocking.py -f Vinardo -g SampleACE2RefLigand.pdb
1601 -r SampleACE2Receptor.pdbqt -i SampleACE2Ligands.sdf
1602 -o SampleACE2LigandsOut.sdf
1603
1604 To run the first example for docking molecules using AD4 forcefield relying on
1605 the presence of affinity maps in the working directory and write out a SD file,
1606 type:
1607
1608 % VinaPerformDocking.py -f AD4 -g SampleACE2RefLigand.pdb
1609 -r SampleACE2Receptor -i SampleACE2Ligands.sdf
1610 -o SampleACE2LigandsOut.sdf
1611
1612 To run the first example for docking molecules using a set of explicit values
1613 for grid dimensions and write out a SD file, type:
1614
1615 % VinaPerformDocking.py -g "41.399, 5.851, 28.082"
1616 --gridSize "25.0, 25.0, 25.0" --gridSpacing 0.375
1617 -r SampleACE2Receptor.pdbqt -i SampleACE2Ligands.sdf
1618 -o SampleACE2LigandsOut.sdf
1619
1620 To run the first example for only scoring molecules already positioned in the
1621 binding pocket and write out a SD file, type:
1622
1623 % VinaPerformDocking.py -m ScoreOnly -g SampleACE2RefLigand.sdf
1624 -r SampleACE2Receptor.pdbqt -i SampleACE2RefLigandWithHs.sdf
1625 -o SampleACE2LigandsOut.sdf
1626
1627 To run the first example for docking molecules to increase the accuracy of
1628 pose predictions and write out a SD file, type:
1629
1630 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1631 -r SampleACE2Receptor.pdbqt --exhaustiveness 24
1632 -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1633
1634 To run the first example in multiprocessing mode on all available CPUs
1635 without loading all data into memory, a single thread for Vina docking, and
1636 write out a SD file, type:
1637
1638 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1639 -r SampleACE2Receptor.pdbqt --mp yes
1640 -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1641
1642 To run the first example in multiprocessing mode on all available CPUs
1643 by loading all data into memory, a single thread for Vina, and write out
1644 a SD file, type:
1645
1646 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1647 -r SampleACE2Receptor.pdbqt --mp yes --mpParams "inputDataMode,
1648 InMemory" -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1649
1650 To run the first example in multiprocessing mode on specific number of CPUs
1651 and chunk size without loading all data into memory along with a specific number
1652 of threads for Vina docking and write out a SD file, type:
1653
1654 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1655 -r SampleACE2Receptor.pdbqt --mp yes --mpParams "inputDataMode,lazy,
1656 numProcesses,4,chunkSize,2" --numThreads 2
1657 -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1658
1659 To run the first example for docking molecules employing a flexible portion
1660 of the receptor corresponding to ARG273 and write out a SD file, type:
1661
1662 % VinaPerformDocking.py -g SampleACE2RefLigand.pdb
1663 -r SampleACE2RigidReceptor.pdbqt
1664 --receptorFlexFile SampleACE2FlexReceptor.pdbqt
1665 -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1666
1667 To run the first example for docking molecules using specified parameters and
1668 write out a SD file, type:
1669
1670 % VinaPerformDocking.py -g "41.399, 5.851, 28.082"
1671 --gridSize "25.0, 25.0, 25.0" --gridSpacing 0.375 --energyComponents
1672 yes --exhaustiveness 32 --forcefield Vina --mode dock --numPoses 2
1673 --numThreads 4 --randomSeed 42 --validateMolecules no
1674 --vinaVerbosity 0 -r SampleACE2Receptor.pdbqt
1675 -i SampleACE2Ligands.sdf -o SampleACE2LigandsOut.sdf
1676
1677 Author:
1678 Manish Sud(msud@san.rr.com)
1679
1680 Acknowledgments:
1681 Diogo Santos-Martins and Stefano Forli
1682
1683 See also:
1684 PyMOLConvertLigandFileFormat.py, PyMOLExtractSelection.py,
1685 PyMOLInfoMacromolecules.py, PyMOLVisualizeMacromolecules.py,
1686 RDKitConvertFileFormat.py, RDKitEnumerateTautomers.py,
1687 RDKitGenerateConformers.py, RDKitPerformMinimization.py,
1688 RDKitPerformConstrainedMinimization.py
1689
1690 Copyright:
1691 Copyright (C) 2026 Manish Sud. All rights reserved.
1692
1693 The functionality available in this script is implemented using AutoDockVina
1694 and Meeko, open source software packages for docking, and RDKit, an open
1695 source toolkit for cheminformatics developed by Greg Landrum.
1696
1697 This file is part of MayaChemTools.
1698
1699 MayaChemTools is free software; you can redistribute it and/or modify it under
1700 the terms of the GNU Lesser General Public License as published by the Free
1701 Software Foundation; either version 3 of the License, or (at your option) any
1702 later version.
1703
1704 """
1705
1706 if __name__ == "__main__":
1707 main()
