1 #!/bin/env python
   2 #
   3 # File: OpenMMPrepareMacromolecule.py
   4 # Author: Manish Sud <msud@san.rr.com>
   5 #
   6 # Copyright (C) 2026 Manish Sud. All rights reserved.
   7 #
   8 # The functionality available in this script is implemented using OpenMM, an
   9 # open source molecuar simulation package.
  10 #
  11 # This file is part of MayaChemTools.
  12 #
  13 # MayaChemTools is free software; you can redistribute it and/or modify it under
  14 # the terms of the GNU Lesser General Public License as published by the Free
  15 # Software Foundation; either version 3 of the License, or (at your option) any
  16 # later version.
  17 #
  18 # MayaChemTools is distributed in the hope that it will be useful, but without
  19 # any warranty; without even the implied warranty of merchantability of fitness
  20 # for a particular purpose.  See the GNU Lesser General Public License for more
  21 # details.
  22 #
  23 # You should have received a copy of the GNU Lesser General Public License
  24 # along with MayaChemTools; if not, see <http://www.gnu.org/licenses/> or
  25 # write to the Free Software Foundation Inc., 59 Temple Place, Suite 330,
  26 # Boston, MA, 02111-1307, USA.
  27 #
  28 
  29 from __future__ import print_function
  30 
  31 import os
  32 import sys
  33 import time
  34 import re
  35 
  36 # OpenMM imports...
  37 try:
  38     import openmm as mm
  39     import pdbfixer
  40     from pdbfixer.pdbfixer import __version__ as pdbfixerversion
  41 except ImportError as ErrMsg:
  42     sys.stderr.write("\nFailed to import OpenMM related module/package: %s\n" % ErrMsg)
  43     sys.stderr.write("Check/update your OpenMM environment and try again.\n\n")
  44     sys.exit(1)
  45 
  46 # MayaChemTools imports...
  47 sys.path.insert(0, os.path.join(os.path.dirname(sys.argv[0]), "..", "lib", "Python"))
  48 try:
  49     from docopt import docopt
  50     import MiscUtil
  51     import OpenMMUtil
  52 except ImportError as ErrMsg:
  53     sys.stderr.write("\nFailed to import MayaChemTools module/package: %s\n" % ErrMsg)
  54     sys.stderr.write("Check/update your MayaChemTools environment and try again.\n\n")
  55     sys.exit(1)
  56 
  57 ScriptName = os.path.basename(sys.argv[0])
  58 Options = {}
  59 OptionsInfo = {}
  60 
  61 
  62 def main():
  63     """Start execution of the script."""
  64 
  65     MiscUtil.PrintInfo(
  66         "\n%s (OpenMM v%s; PDBFixerVersion v%s; MayaChemTools v%s; %s): Starting...\n"
  67         % (
  68             ScriptName,
  69             mm.Platform.getOpenMMVersion(),
  70             pdbfixerversion,
  71             MiscUtil.GetMayaChemToolsVersion(),
  72             time.asctime(),
  73         )
  74     )
  75 
  76     (WallClockTime, ProcessorTime) = MiscUtil.GetWallClockAndProcessorTime()
  77 
  78     # Retrieve command line arguments and options...
  79     RetrieveOptions()
  80 
  81     # Process and validate command line arguments and options...
  82     ProcessOptions()
  83 
  84     # Perform actions required by the script...
  85     PrepareMacromolecule()
  86 
  87     MiscUtil.PrintInfo("\n%s: Done...\n" % ScriptName)
  88     MiscUtil.PrintInfo("Total time: %s" % MiscUtil.GetFormattedElapsedTime(WallClockTime, ProcessorTime))
  89 
  90 
  91 def PrepareMacromolecule():
  92     """Prepare a macromolecule for simulation and write it out."""
  93 
  94     # Read macromolecule...
  95     PDBFixerObjectHandle = ReadMacromolecule()
  96 
  97     # Identify missing residues...
  98     IdentifyMissingResidues(PDBFixerObjectHandle)
  99 
 100     # Identify and replace non-standard residues...
 101     IdentifyAndReplaceNonStandardResidues(PDBFixerObjectHandle)
 102 
 103     # Identify and delete heterogen residues...
 104     IdentifyAndDeleteHeterogenResidues(PDBFixerObjectHandle)
 105 
 106     # Identify and add missing atoms...
 107     IdentifyAndAddMissingAtoms(PDBFixerObjectHandle)
 108 
 109     if OptionsInfo["WaterBox"]:
 110         AddWaterBox(PDBFixerObjectHandle)
 111     elif OptionsInfo["Membrane"]:
 112         AddLipidMembrane(PDBFixerObjectHandle)
 113     else:
 114         MiscUtil.PrintInfo("\nSkipping addtion of water box or membrane...")
 115 
 116     # Write macromolecule...
 117     WriteMacromolecule(PDBFixerObjectHandle)
 118 
 119 
 120 def ReadMacromolecule():
 121     """Read macromolecule."""
 122 
 123     Infile = OptionsInfo["Infile"]
 124     MiscUtil.PrintInfo("\nProcessing file %s..." % Infile)
 125 
 126     PDBFixerObjectHandle = pdbfixer.pdbfixer.PDBFixer(filename=Infile)
 127 
 128     return PDBFixerObjectHandle
 129 
 130 
 131 def WriteMacromolecule(PDBFixerObjectHandle):
 132     """Write macromolecule."""
 133 
 134     MiscUtil.PrintInfo("\nGenerating output file %s..." % OptionsInfo["Outfile"])
 135     OpenMMUtil.WritePDBFile(
 136         OptionsInfo["Outfile"], PDBFixerObjectHandle.topology, PDBFixerObjectHandle.positions, KeepIDs=True
 137     )
 138 
 139 
 140 def IdentifyMissingResidues(PDBFixerObjectHandle):
 141     """Identify missing residues based on PDB REMARK records."""
 142 
 143     MiscUtil.PrintInfo("\nIdentifying missing residues...")
 144 
 145     PDBFixerObjectHandle.missingResidues = {}
 146     PDBFixerObjectHandle.findMissingResidues()
 147 
 148     ListMissingResidues(PDBFixerObjectHandle)
 149 
 150 
 151 def ListMissingResidues(PDBFixerObjectHandle):
 152     """List missing residues."""
 153 
 154     if len(PDBFixerObjectHandle.missingResidues) == 0:
 155         MiscUtil.PrintInfo("Number of missing residues: 0")
 156         return
 157 
 158     Chains = list(PDBFixerObjectHandle.topology.chains())
 159 
 160     MissingResiduesCount = 0
 161     MissingResiduesInfo = []
 162     for Index, Key in enumerate(sorted(PDBFixerObjectHandle.missingResidues)):
 163         ChainIndex = Key[0]
 164         # ResidueIndex corresponds to the residue after the insertion of missing residues...`
 165         ResidueIndex = Key[1]
 166 
 167         MissingResidues = PDBFixerObjectHandle.missingResidues[Key]
 168         MissingResiduesCount += len(MissingResidues)
 169 
 170         Chain = Chains[ChainIndex]
 171         ChainResidues = list(Chain.residues())
 172 
 173         if ResidueIndex < len(ChainResidues):
 174             ResidueNumOffset = int(ChainResidues[ResidueIndex].id) - len(MissingResidues) - 1
 175         else:
 176             ResidueNumOffset = int(ChainResidues[-1].id)
 177 
 178         StartResNum = ResidueNumOffset + 1
 179         EndResNum = ResidueNumOffset + len(MissingResidues)
 180 
 181         MissingResiduesInfo.append(
 182             "Chain: %s; StartResNum-EndResNum: %s-%s; Count: %s; Residues: %s"
 183             % (Chain.id, StartResNum, EndResNum, len(MissingResidues), ", ".join(MissingResidues))
 184         )
 185 
 186     MiscUtil.PrintInfo("Total number of missing residues: %s" % (MissingResiduesCount))
 187     if OptionsInfo["ListDetails"]:
 188         MiscUtil.PrintInfo("%s" % ("\n".join(MissingResiduesInfo)))
 189 
 190 
 191 def IdentifyAndReplaceNonStandardResidues(PDBFixerObjectHandle):
 192     """Identify and replace non-standard residues."""
 193 
 194     MiscUtil.PrintInfo("\nIdentifying non-standard residues...")
 195 
 196     PDBFixerObjectHandle.findNonstandardResidues()
 197     NonStandardResidues = PDBFixerObjectHandle.nonstandardResidues
 198     NumNonStandardResiues = len(NonStandardResidues)
 199 
 200     MiscUtil.PrintInfo("Total number of non-standard residues: %s" % NumNonStandardResiues)
 201 
 202     if NumNonStandardResiues == 0:
 203         return
 204 
 205     if OptionsInfo["ListDetails"]:
 206         ListNonStandardResidues(PDBFixerObjectHandle)
 207 
 208     if not OptionsInfo["ReplaceNonStandardResidues"]:
 209         MiscUtil.PrintInfo("Skipping replacement of non-standard residues...")
 210         return
 211 
 212     MiscUtil.PrintInfo("Replacing non-standard residues...")
 213     PDBFixerObjectHandle.replaceNonstandardResidues()
 214 
 215 
 216 def ListNonStandardResidues(PDBFixerObjectHandle):
 217     """List non-standard residues."""
 218 
 219     NonStandardResidues = PDBFixerObjectHandle.nonstandardResidues
 220 
 221     if len(NonStandardResidues) == 0:
 222         return
 223 
 224     NonStandardResiduesMappingInfo = []
 225     for Values in NonStandardResidues:
 226         NonStandardRes = Values[0]
 227         StandardResName = Values[1]
 228 
 229         MappingInfo = "<%s %s %s>: %s" % (
 230             NonStandardRes.id,
 231             NonStandardRes.name,
 232             NonStandardRes.chain.id,
 233             StandardResName,
 234         )
 235         NonStandardResiduesMappingInfo.append(MappingInfo)
 236 
 237     if len(NonStandardResiduesMappingInfo):
 238         MiscUtil.PrintInfo("Non-standard residues mapping: %s\n" % ",".join(NonStandardResiduesMappingInfo))
 239 
 240 
 241 def IdentifyAndDeleteHeterogenResidues(PDBFixerObjectHandle):
 242     """Identify and delete heterogen residues."""
 243 
 244     DeleteHeterogens = OptionsInfo["DeleteHeterogens"]
 245     if DeleteHeterogens is None:
 246         MiscUtil.PrintInfo("\nSkipping deletion of any heterogen residues...")
 247         return
 248 
 249     if re.match("^All$", DeleteHeterogens, re.I):
 250         MiscUtil.PrintInfo("\nDeleting all heterogen residues including water...")
 251         PDBFixerObjectHandle.removeHeterogens(keepWater=False)
 252     elif re.match("^AllExceptWater$", DeleteHeterogens, re.I):
 253         MiscUtil.PrintInfo("\nDeleting all heterogen residues except water...")
 254         PDBFixerObjectHandle.removeHeterogens(keepWater=True)
 255     elif re.match("^WaterOnly$", DeleteHeterogens, re.I):
 256         MiscUtil.PrintInfo("\nDeleting water only heterogen residues...")
 257         DeleteWater(PDBFixerObjectHandle)
 258     else:
 259         MiscUtil.PrintError(
 260             'The value specified, %s, for option "--deleteHeterogens" is not valid. Supported values: All AllExceptWater WaterOnly None'
 261             % DeleteHeterogens
 262         )
 263 
 264 
 265 def DeleteWater(PDBFixerObjectHandle):
 266     """Delete water."""
 267 
 268     ModellerObjectHandle = mm.app.Modeller(PDBFixerObjectHandle.topology, PDBFixerObjectHandle.positions)
 269     ModellerObjectHandle.deleteWater()
 270 
 271     PDBFixerObjectHandle.topology = ModellerObjectHandle.topology
 272     PDBFixerObjectHandle.positions = ModellerObjectHandle.positions
 273 
 274 
 275 def IdentifyAndAddMissingAtoms(PDBFixerObjectHandle):
 276     """Identify and missing atoms along with already identified missing residues."""
 277 
 278     MiscUtil.PrintInfo("\nIdentifying missing atoms...")
 279     PDBFixerObjectHandle.findMissingAtoms()
 280 
 281     ListMissingAtoms(PDBFixerObjectHandle)
 282 
 283     if OptionsInfo["AddHeavyAtoms"] and OptionsInfo["AddResidues"]:
 284         Msg = "Adding missing heavy atoms along with missing residues..."
 285     elif OptionsInfo["AddHeavyAtoms"]:
 286         Msg = "Adding missing heavy atoms..."
 287     elif OptionsInfo["AddResidues"]:
 288         Msg = "Adding missing residues..."
 289     else:
 290         Msg = "Skipping addition of any heavy atoms along with any missing residues..."
 291 
 292     if not OptionsInfo["AddHeavyAtoms"]:
 293         PDBFixerObjectHandle.missingAtoms = {}
 294         PDBFixerObjectHandle.missingTerminals = {}
 295 
 296     MiscUtil.PrintInfo("\n%s" % Msg)
 297     if OptionsInfo["AddHeavyAtoms"] or OptionsInfo["AddResidues"]:
 298         PDBFixerObjectHandle.addMissingAtoms()
 299 
 300     if OptionsInfo["AddHydrogens"]:
 301         MiscUtil.PrintInfo("Adding missing hydrogens at pH %s..." % OptionsInfo["AddHydrogensAtpH"])
 302         PDBFixerObjectHandle.addMissingHydrogens(pH=OptionsInfo["AddHydrogensAtpH"], forcefield=None)
 303     else:
 304         MiscUtil.PrintInfo("Skipping addition of any missing hydrogens...")
 305 
 306 
 307 def ListMissingAtoms(PDBFixerObjectHandle):
 308     """List missing atoms."""
 309 
 310     if len(PDBFixerObjectHandle.missingAtoms) == 0 and len(PDBFixerObjectHandle.missingTerminals) == 0:
 311         MiscUtil.PrintInfo("Total number of missing atoms: 0")
 312         return
 313 
 314     ListMissingAtomsInfo(PDBFixerObjectHandle.missingAtoms, TerminalAtoms=False)
 315     ListMissingAtomsInfo(PDBFixerObjectHandle.missingTerminals, TerminalAtoms=True)
 316 
 317 
 318 def ListMissingAtomsInfo(MissingAtoms, TerminalAtoms):
 319     """Setup missing atoms information."""
 320 
 321     MissingAtomsInfo = []
 322     MissingAtomsCount = 0
 323     MissingAtomsResiduesCount = 0
 324     for Residue, Atoms in MissingAtoms.items():
 325         MissingAtomsResiduesCount += 1
 326         MissingAtomsCount += len(Atoms)
 327 
 328         if TerminalAtoms:
 329             AtomsNames = [AtomName for AtomName in Atoms]
 330         else:
 331             AtomsNames = [Atom.name for Atom in Atoms]
 332         AtomsInfo = "<%s %s %s>: %s" % (Residue.id, Residue.name, Residue.chain.id, ", ".join(AtomsNames))
 333         MissingAtomsInfo.append(AtomsInfo)
 334 
 335     AtomsLabel = "terminal atoms" if TerminalAtoms else "atoms"
 336     if MissingAtomsCount == 0:
 337         MiscUtil.PrintInfo("Total number of missing %s: %s" % (MissingAtomsCount, AtomsLabel))
 338     else:
 339         ResiduesLabel = "residue" if MissingAtomsResiduesCount == 1 else "residues"
 340         MiscUtil.PrintInfo(
 341             "Total number of missing %s across %s %s: %s"
 342             % (AtomsLabel, MissingAtomsResiduesCount, ResiduesLabel, MissingAtomsCount)
 343         )
 344         if OptionsInfo["ListDetails"]:
 345             MiscUtil.PrintInfo("Missing %s across residues: %s" % (AtomsLabel, "; ".join(MissingAtomsInfo)))
 346 
 347 
 348 def AddWaterBox(PDBFixerObjectHandle):
 349     """Add water box."""
 350 
 351     MiscUtil.PrintInfo("\nAdding water box...")
 352 
 353     WaterBoxParams = OptionsInfo["WaterBoxParams"]
 354 
 355     Size, Padding, Shape = [None] * 3
 356     if WaterBoxParams["ModeSize"]:
 357         SizeList = WaterBoxParams["SizeList"]
 358         Size = mm.Vec3(SizeList[0], SizeList[1], SizeList[2]) * mm.unit.nanometer
 359     elif WaterBoxParams["ModePadding"]:
 360         Padding = WaterBoxParams["Padding"] * mm.unit.nanometer
 361         Shape = WaterBoxParams["Shape"]
 362     else:
 363         MiscUtil.PrintError(
 364             'The parameter value, %s, specified for parameter name, mode, using "--waterBoxParams" option is not a valid value. Supported values: Size Padding \n'
 365             % (WaterBoxParams["Mode"])
 366         )
 367 
 368     IonicStrength = OptionsInfo["IonicStrength"] * mm.unit.molar
 369     PDBFixerObjectHandle.addSolvent(
 370         boxSize=Size,
 371         padding=Padding,
 372         boxShape=Shape,
 373         positiveIon=OptionsInfo["IonPositive"],
 374         negativeIon=OptionsInfo["IonNegative"],
 375         ionicStrength=IonicStrength,
 376     )
 377 
 378 
 379 def AddLipidMembrane(PDBFixerObjectHandle):
 380     """Add lipid membrane along with water."""
 381 
 382     MiscUtil.PrintInfo("\nAdding membrane along with water...")
 383 
 384     MembraneParams = OptionsInfo["MembraneParams"]
 385 
 386     LipidType = MembraneParams["LipidType"]
 387     MembraneCenterZ = MembraneParams["MembraneCenterZ"] * mm.unit.nanometer
 388     Padding = MembraneParams["Padding"] * mm.unit.nanometer
 389 
 390     IonicStrength = OptionsInfo["IonicStrength"] * mm.unit.molar
 391 
 392     PDBFixerObjectHandle.addMembrane(
 393         lipidType=LipidType,
 394         membraneCenterZ=MembraneCenterZ,
 395         minimumPadding=Padding,
 396         positiveIon=OptionsInfo["IonPositive"],
 397         negativeIon=OptionsInfo["IonNegative"],
 398         ionicStrength=IonicStrength,
 399     )
 400 
 401 
 402 def ProcessMembraneParamsOption():
 403     """Process option for membrane parameters."""
 404 
 405     ParamsOptionName = "--membraneParams"
 406     ParamsOptionValue = Options[ParamsOptionName]
 407     ParamsDefaultInfo = {"LipidType": ["str", "POPC"], "MembraneCenterZ": ["float", 0.0], "Padding": ["float", 1.0]}
 408     MembraneParams = MiscUtil.ProcessOptionNameValuePairParameters(
 409         ParamsOptionName, ParamsOptionValue, ParamsDefaultInfo
 410     )
 411 
 412     for ParamName in ["LipidType", "MembraneCenterZ", "Padding"]:
 413         ParamValue = MembraneParams[ParamName]
 414         if ParamName == "LipidType":
 415             LipidTypes = ["POPC", "POPE", "DLPC", "DLPE", "DMPC", "DOPC", "DPPC"]
 416             if ParamValue not in LipidTypes:
 417                 MiscUtil.PrintError(
 418                     'The parameter value, %s, specified for parameter name, %s, using "%s" option is not a valid value. Supported values: %s \n'
 419                     % (ParamValue, ParamName, ParamsOptionName, LipidTypes)
 420                 )
 421         elif ParamName == "Padding":
 422             if ParamValue <= 0:
 423                 MiscUtil.PrintError(
 424                     'The parameter value, %s, specified for parameter name, %s, using "%s" option is not a valid value. Supported values: > 0\n'
 425                     % (ParamValue, ParamName, ParamsOptionName)
 426                 )
 427 
 428     OptionsInfo["MembraneParams"] = MembraneParams
 429 
 430 
 431 def ProcessOptions():
 432     """Process and validate command line arguments and options."""
 433 
 434     MiscUtil.PrintInfo("Processing options...")
 435 
 436     # Validate options...
 437     ValidateOptions()
 438 
 439     OptionsInfo["Infile"] = Options["--infile"]
 440     FileDir, FileName, FileExt = MiscUtil.ParseFileName(OptionsInfo["Infile"])
 441     OptionsInfo["InfileRoot"] = FileName
 442 
 443     OptionsInfo["Outfile"] = Options["--outfile"]
 444 
 445     OptionsInfo["AddHeavyAtoms"] = True if re.match("^yes$", Options["--addHeavyAtoms"], re.I) else False
 446 
 447     OptionsInfo["AddHydrogens"] = True if re.match("^yes$", Options["--addHydrogens"], re.I) else False
 448     OptionsInfo["AddHydrogensAtpH"] = float(Options["--addHydrogensAtpH"])
 449 
 450     OptionsInfo["AddResidues"] = True if re.match("^yes$", Options["--addResidues"], re.I) else False
 451 
 452     DeleteHeterogens = Options["--deleteHeterogens"]
 453     if re.match("^All$", DeleteHeterogens, re.I):
 454         DeleteHeterogens = "All"
 455     elif re.match("^AllExceptWater$", DeleteHeterogens, re.I):
 456         DeleteHeterogens = "AllExceptWater"
 457     elif re.match("^WaterOnly$", DeleteHeterogens, re.I):
 458         DeleteHeterogens = "WaterOnly"
 459     elif re.match("^None$", DeleteHeterogens, re.I):
 460         DeleteHeterogens = None
 461     elif re.match("^auto$", DeleteHeterogens, re.I):
 462         DeleteHeterogens = None
 463         if re.match("^yes$", Options["--membrane"], re.I) or re.match("^yes$", Options["--waterBox"], re.I):
 464             DeleteHeterogens = "WaterOnly"
 465     else:
 466         MiscUtil.PrintError(
 467             'The value specified, %s, for option "--deleteHeterogens" is not valid. Supported values: All AllExceptWater WaterOnly None'
 468             % DeleteHeterogens
 469         )
 470     OptionsInfo["DeleteHeterogens"] = DeleteHeterogens
 471 
 472     OptionsInfo["IonPositive"] = Options["--ionPositive"]
 473     OptionsInfo["IonNegative"] = Options["--ionNegative"]
 474     OptionsInfo["IonicStrength"] = float(Options["--ionicStrength"])
 475 
 476     OptionsInfo["ListDetails"] = True if re.match("^yes$", Options["--listDetails"], re.I) else False
 477 
 478     OptionsInfo["Membrane"] = True if re.match("^yes$", Options["--membrane"], re.I) else False
 479     ProcessMembraneParamsOption()
 480 
 481     OptionsInfo["ReplaceNonStandardResidues"] = (
 482         True if re.match("^yes$", Options["--replaceNonStandardResidues"], re.I) else False
 483     )
 484 
 485     OptionsInfo["WaterBox"] = True if re.match("^yes$", Options["--waterBox"], re.I) else False
 486     OptionsInfo["WaterBoxParams"] = OpenMMUtil.ProcessOptionOpenMMWaterBoxParameters(
 487         "--waterBoxParams", Options["--waterBoxParams"]
 488     )
 489 
 490     OptionsInfo["Overwrite"] = Options["--overwrite"]
 491 
 492 
 493 def RetrieveOptions():
 494     """Retrieve command line arguments and options."""
 495 
 496     # Get options...
 497     global Options
 498     Options = docopt(_docoptUsage_)
 499 
 500     # Set current working directory to the specified directory...
 501     WorkingDir = Options["--workingdir"]
 502     if WorkingDir:
 503         os.chdir(WorkingDir)
 504 
 505     # Handle examples option...
 506     if "--examples" in Options and Options["--examples"]:
 507         MiscUtil.PrintInfo(MiscUtil.GetExamplesTextFromDocOptText(_docoptUsage_))
 508         sys.exit(0)
 509 
 510 
 511 def ValidateOptions():
 512     """Validate option values."""
 513 
 514     MiscUtil.ValidateOptionFilePath("-i, --infile", Options["--infile"])
 515     MiscUtil.ValidateOptionFileExt("-i, --infile", Options["--infile"], "pdb cif")
 516 
 517     MiscUtil.ValidateOptionFileExt("-o, --outfile", Options["--outfile"], "pdb cif")
 518     MiscUtil.ValidateOptionsOutputFileOverwrite(
 519         "-o, --outfile", Options["--outfile"], "--overwrite", Options["--overwrite"]
 520     )
 521     MiscUtil.ValidateOptionsDistinctFileNames(
 522         "-i, --infile", Options["--infile"], "-o, --outfile", Options["--outfile"]
 523     )
 524 
 525     MiscUtil.ValidateOptionTextValue("--addHeavyAtoms", Options["--addHeavyAtoms"], "yes no")
 526 
 527     MiscUtil.ValidateOptionTextValue("--addHydrogens", Options["--addHydrogens"], "yes no")
 528     MiscUtil.ValidateOptionFloatValue("--addHydrogensAtpH", Options["--addHydrogensAtpH"], {">=": 0})
 529 
 530     MiscUtil.ValidateOptionTextValue("--addResidues", Options["--addResidues"], "yes no")
 531 
 532     MiscUtil.ValidateOptionTextValue(
 533         "--deleteHeterogens", Options["--deleteHeterogens"], "All AllExceptWater WaterOnly None auto"
 534     )
 535     DeleteHeterogens = Options["--deleteHeterogens"]
 536     if re.match("^(AllExceptWater|None)$", DeleteHeterogens, re.I):
 537         if re.match("^yes$", Options["--membrane"], re.I):
 538             MiscUtil.PrintError(
 539                 'The value specified, %s, for option "--deleteHeterogens" is not valid during "yes" value of option "--membrane". '
 540                 % DeleteHeterogens
 541             )
 542         elif re.match("^yes$", Options["--waterBox"], re.I):
 543             MiscUtil.PrintError(
 544                 'The value specified, %s, for option "--deleteHeterogens" is not valid during "yes" value of option "--waterBox". '
 545                 % DeleteHeterogens
 546             )
 547 
 548     ValidValues = "Li+ Na+ K+ Rb+ Cs+"
 549     EscapedValidValuesPattern = r"Li\+|Na\+|K\+|Rb\+|Cs\+"
 550     Value = Options["--ionPositive"]
 551     if not re.match("^(%s)$" % EscapedValidValuesPattern, Value):
 552         MiscUtil.PrintError(
 553             'The value specified, %s, for option "-ionPositive" is not valid: Supported value(s): %s'
 554             % (Value, ValidValues)
 555         )
 556 
 557     ValidValues = "F- Cl- Br- I-"
 558     ValidValuesPattern = "F-|Cl-|Br-|I-"
 559     Value = Options["--ionNegative"]
 560     if not re.match("^(%s)$" % ValidValuesPattern, Value):
 561         MiscUtil.PrintError(
 562             'The value specified, %s, for option "-ionNegative" is not valid: Supported value(s): %s'
 563             % (Value, ValidValues)
 564         )
 565 
 566     MiscUtil.ValidateOptionFloatValue("--ionicStrength", Options["--ionicStrength"], {">=": 0})
 567 
 568     MiscUtil.ValidateOptionTextValue("--listDetails", Options["--listDetails"], "yes no")
 569     MiscUtil.ValidateOptionTextValue("--replaceNonStandardResidues", Options["--replaceNonStandardResidues"], "yes no")
 570 
 571     MiscUtil.ValidateOptionTextValue("--membrane", Options["--membrane"], "yes no")
 572     MiscUtil.ValidateOptionTextValue("--waterBox", Options["--waterBox"], "yes no")
 573 
 574     if re.match("^yes$", Options["--membrane"], re.I) and re.match("^yes$", Options["--waterBox"], re.I):
 575         MiscUtil.PrintError(
 576             'You\'ve specified, "Yes" for both "--membrane" and "--waterBox" options. It\'s not allowed. You must choose between adding a water box or a membrane. The water is automatically added during the construction of the membrane.'
 577         )
 578 
 579 
 580 # Setup a usage string for docopt...
 581 _docoptUsage_ = """
 582 OpenMMPrepareMacromolecule.py - Prepare a macromolecule for simulation
 583 
 584 Usage:
 585     OpenMMPrepareMacromolecule.py [--addHeavyAtoms <yes or no>] [--addHydrogens <yes or no>] [--addHydrogensAtpH <number>]
 586                                   [--addResidues <yes or no>] [--deleteHeterogens <All, AllExceptWater, WaterOnly, None>] [--ionPositive <text>]
 587                                   [--ionNegative <text>] [--ionicStrength <number>] [--listDetails <yes or no> ] [--membrane <yes or no>]
 588                                   [--membraneParams <Name,Value,..>] [--overwrite] [--replaceNonStandardResidues <yes or no>] [--waterBox <yes or no>]
 589                                   [--waterBoxParams <Name,Value,..>] [-w <dir>] -i <infile> -o <outfile>
 590     OpenMMPrepareMacromolecule.py -h | --help | -e | --examples
 591 
 592 Description:
 593     Prepare a macromolecute in an input file for molecular simulation and
 594     and write it out to an output file. The macromolecule is prepared by
 595     automatically performing the following actions:
 596 
 597         . Identify and replace non-standard residues
 598         . Add missing residues
 599         . Add missing heavy atoms
 600         . Add missing hydrogens
 601 
 602     You may optionally remove heterogens, add a water box, and add a lipid
 603     membrane.
 604 
 605     The supported input file format are:  PDB (.pdb) and CIF (.cif)
 606 
 607     The supported output file formats are:  PDB (.pdb) and CIF (.cif)
 608 
 609 Options:
 610     --addHeavyAtoms <yes or no>  [default: yes]
 611         Add missing non-hydrogen atoms based on the templates. The terminal atoms
 612         are also added.
 613     --addHydrogens <yes or no>  [default: yes]
 614         Add missing hydrogens at pH value specified using '-addHydrogensAtpH'
 615         option. The missing hydrogens are added based on the templates and pKa
 616         calculations are performed.
 617     --addHydrogensAtpH <number>  [default: 7.0]
 618         pH value to use for adding missing hydrogens.
 619     --addResidues <yes or no>  [default: yes]
 620         Add missing residues unidentified based on the PDB records.
 621     --deleteHeterogens <All, AllExceptWater, WaterOnly, None>  [default: auto]
 622         Delete heterogens corresponding to  non-standard names of amino acids, dna,
 623         and rna along with any ligand names. 'N' and 'UNK' also consider standard
 624         residues. Default value: WaterOnly during addition of WaterBox or Membrane;
 625         Otherwise, None.
 626         
 627         The 'AllExceptWater' or 'None' values are not allowed during the addition of
 628         a water box or membrane. The waters must be deleted as they are explicitly
 629         added during the construction of a water box and membrane.
 630     -e, --examples
 631         Print examples.
 632     -h, --help
 633         Print this help message.
 634     -i, --infile <infile>
 635         Input file name.
 636     --ionPositive <text>  [default: Na+]
 637         Type of positive ion to add during the addition of a water box or membrane.
 638         Possible values: Li+, Na+, K+, Rb+, or Cs+.
 639     --ionNegative <text>  [default: Cl-]
 640         Type of negative ion to add during the addition of a water box or membrane.
 641         Possible values: Cl-, Br-, F-, or I-.
 642     --ionicStrength <number>  [default: 0.0]
 643         Total concentration of both positive and negative ions to add excluding
 644         the ions added to neutralize the system during the addition of a water box
 645         or a membrane.
 646     -l, --listDetails <yes or no>  [default: no]
 647         List details about missing and non-standard residues along with residues
 648         containing missing atoms.
 649     --membrane <yes or no>  [default: no]
 650         Add lipid membrane along with a water box. The script replies on OpenMM
 651         method modeller.addMebrane() to perform the task. The membrane is added
 652         in the XY plane. The existing macromolecule must be oriented and positioned
 653         correctly.
 654         
 655         A word to the wise: You may want to start with a model from the Orientations
 656         of Proteins in Membranes (OPM) database at http://opm.phar.umich.edu.
 657         
 658         The size of the membrane and water box are determined by the value of
 659         'padding' parameter specified using '--membraneParams' option. All atoms
 660         in macromolecule are guaranteed to be at least this far from any edge of
 661         the periodic box.
 662     --membraneParams <Name,Value,..>  [default: auto]
 663         A comma delimited list of parameter name and value pairs for adding
 664         a lipid membrane and water.
 665         
 666         The supported parameter names along with their default values are
 667         are shown below:
 668             
 669             lipidType, POPC  [ Possible values: POPC, POPE, DLPC, DLPE,  DMPC, 
 670                 DOPC or DPPC ]
 671             membraneCenterZ, 0.0
 672             padding, 1.0
 673             
 674         A brief description of parameters is provided below:
 675             
 676             lipidType: Type of lipid to use for constructing the membrane.
 677             membraneCenterZ: Position along the Z axis of the center of
 678                 the membrane in nanomertes.
 679             padding: Minimum padding distance to use in nanomertes. It's used
 680                 to determine the size of the membrane and water box. All atoms
 681                 in the macromolecule are guaranteed to be at least this far from
 682                 any edge of the periodic box.
 683             
 684     -o, --outfile <outfile>
 685         Output file name.
 686     --overwrite
 687         Overwrite existing files.
 688     --replaceNonStandardResidues <yes or no>  [default: yes]
 689         Replace non-standard residue names by standard residue names based on the
 690         list of non-standard residues available in pdbfixer.
 691     --waterBox <yes or no>  [default: no]
 692         Add water box.
 693     --waterBoxParams <Name,Value,..>  [default: auto]
 694         A comma delimited list of parameter name and value pairs for adding
 695         a water box.
 696         
 697         The supported parameter names along with their default values are
 698         are shown below:
 699             
 700             mode, Padding  [ Possible values: Size or Padding ]
 701             size, None  [ Possible value: xsize ysize zsize ]
 702             padding, 1.0
 703             shape, cube  [ Possible values: cube, dodecahedron, or octahedron ]
 704             
 705         A brief description of parameters is provided below:
 706             
 707             mode: Specify the size of the waterbox explicitly or calculate it
 708                 automatically for a macromolecule along with adding padding
 709                 around ther macromolecule.
 710             size: A space delimited triplet of values corresponding to water
 711                 size in nanometers. It must be specified during 'Size' value of
 712                 'mode' parameter.
 713             padding: Padding around a macromolecule in nanometers for filling
 714                 box with water. It must be specified during 'Padding' value of
 715                 'mode' parameter.
 716             
 717     -w, --workingdir <dir>
 718         Location of working directory which defaults to the current directory.
 719 
 720 Examples:
 721     To prepare a macromolecule in a PDB file by replacing non-standard residues,
 722     adding missing residues, adding missing heavy atoms and missing hydrogens,
 723     and generate a PDB file, type:
 724 
 725         % OpenMMPrepareMacromolecule.py -i Sample11.pdb -o Sample11Out.pdb
 726 
 727     To run the first example for listing additional details about missing atoms and
 728     residues, and generate a PDB file, type:
 729 
 730         % OpenMMPrepareMacromolecule.py --listDetails yes -i Sample11.pdb
 731           -o Sample11Out.pdb
 732 
 733     To run the first example for deleting all heterogens including water along
 734     with performing all default actions, and generate a PDB file, type:
 735 
 736         % OpenMMPrepareMacromolecule.py --deleteHeterogens All -i Sample11.pdb
 737           -o Sample11Out.pdb
 738 
 739     To run the first example for deleting water only heterogens along with performing
 740     all default actions, and generate a PDB file, type:
 741 
 742         % OpenMMPrepareMacromolecule.py --deleteHeterogens WaterOnly
 743           -i Sample11.pdb -o Sample11Out.pdb --ov
 744 
 745     To run the first example for adding a water box by automatically calculating its
 746     size, along with performing all default actions, and generate a PDB file, type:
 747 
 748         % OpenMMPrepareMacromolecule.py --waterBox yes -i Sample11.pdb
 749           -o Sample11Out.pdb
 750 
 751     To run the previous example by explcitly specifying various water box parameters,
 752     and generate a PDB file, type:
 753 
 754         % OpenMMPrepareMacromolecule.py --waterBox yes
 755           --waterBoxParams "mode,Padding, padding, 1.0, shape, cube"
 756           -i Sample11.pdb -o Sample11Out.pdb
 757 
 758     To run the first example for adding a water box of a specific size along with
 759     performing all default actions, and generate a PDB file, type:
 760 
 761         % OpenMMPrepareMacromolecule.py --waterBox yes
 762           --waterBoxParams "mode,Size, size, 7.635 7.077 7.447"
 763           -i Sample11.pdb -o Sample11Out.pdb
 764 
 765     To add a lipid membrane around a correctly oriented and positioned macromolecule,
 766     along with performing all default actions,  and generate a PDB file, type:
 767 
 768         % OpenMMPrepareMacromolecule.py --membrane yes
 769           -i Sample12.pdb -o Sample12Out.pdb
 770 
 771     To add a lipid membrane around a correctly oriented and positioned macromolecule,
 772     deleting all heterogens along with performing all default actions,  and generate a
 773     PDB file, type:
 774 
 775         % OpenMMPrepareMacromolecule.py --membrane yes --deleteHeterogens All 
 776           -i Sample12.pdb -o Sample12Out.pdb
 777 
 778     To run the previous example by explcitly specifying various membrane parameters,
 779     and generate a PDB file, type:
 780 
 781         % OpenMMPrepareMacromolecule.py --membrane yes
 782           --membraneParams "lipidType, POPC, membraneCenterZ, 0.0, padding, 1.0"
 783           -i Sample12.pdb -o Sample12Out.pdb
 784 
 785 Author:
 786     Manish Sud(msud@san.rr.com)
 787 
 788 See also:
 789     InfoPDBFiles.pl, ExtractFromPDBFiles.pl, PyMOLExtractSelection.py,
 790     PyMOLInfoMacromolecules.py, PyMOLSplitChainsAndLigands.py
 791 
 792 Copyright:
 793     Copyright (C) 2026 Manish Sud. All rights reserved.
 794 
 795     The functionality available in this script is implemented using OpenMM, an
 796     open source molecuar simulation package.
 797 
 798     This file is part of MayaChemTools.
 799 
 800     MayaChemTools is free software; you can redistribute it and/or modify it under
 801     the terms of the GNU Lesser General Public License as published by the Free
 802     Software Foundation; either version 3 of the License, or (at your option) any
 803     later version.
 804 
 805 """
 806 
 807 if __name__ == "__main__":
 808     main()