import os import re import shutil import tempfile import commands import math import pyopenmm import simtk.unit as units import openeye.oechem import openeye.oeomega import openeye.oeiupac import openeye.oequacpac def readMolecule(filename, normalize = False): """ Read in a molecule from a file (such as .mol2). ARGUMENTS filename (string) - the name of the file containing a molecule, in a format that OpenEye autodetects (such as .mol2) OPTIONAL ARGUMENTS normalize (boolean) - if True, molecule is normalized (renamed, aromaticity, protonated) after reading (default: False) RETURNS molecule (OEMol) - OEMol representation of molecule EXAMPLES >>> # read a mol2 file >>> molecule = readMolecule('phenol.mol2') >>> # works with any type of file that OpenEye autodetects >>> molecule = readMolecule('phenol.sdf') """ # Open input stream. istream = openeye.oechem.oemolistream() istream.open(filename) # Create molecule. molecule = openeye.oechem.OEMol() # Read the molecule. openeye.oechem.OEReadMolecule(istream, molecule) # Close the stream. istream.close() # Normalize if desired. if normalize: normalizeMolecule(molecule) return molecule def writeMolecule(molecule, filename, substructure_name='MOL', preserve_atomtypes=False): """ Write a molecule to a file in any format OpenEye autodetects from filename (such as .mol2). WARNING: The molecule will be standardized before writing by the high-level OEWriteMolecule function. OEWriteConstMolecule is used, to avoid changing the molecule you pass in. ARGUMENTS molecule (OEMol) - the molecule to be written filename (string) - the file to write the molecule to (type autodetected from filename) OPTIONAL ARGUMENTS substructure_name (String) - if a mol2 file is written, this is used for the substructure name (default: 'MOL') preserve_atomtypes (bool) - if True, a mol2 file will be written with atom types preserved RETURNS None NOTES Multiple conformers are written. EXAMPLES >>> # create a molecule >>> molecule = createMoleculeFromIUPAC('phenol') >>> # write it as a mol2 file >>> writeMolecule(molecule, 'phenol.mol2') """ # Open output stream. ostream = openeye.oechem.oemolostream() ostream.open(filename) # Define internal function for writing multiple conformers to an output stream. def write_all_conformers(ostream, molecule): # write all conformers of each molecule for conformer in molecule.GetConfs(): if preserve_atomtypes: openeye.oechem.OEWriteMol2File(ostream, conformer) else: openeye.oechem.OEWriteConstMolecule(ostream, conformer) return # If 'molecule' is actually a list of molecules, write them all. if type(molecule) == type(list()): for individual_molecule in molecule: write_all_conformers(ostream, individual_molecule) else: write_all_conformers(ostream, molecule) # Close the stream. ostream.close() # Replace substructure name if mol2 file. suffix = os.path.splitext(filename)[-1] if (suffix == '.mol2' and substructure_name != None): modifySubstructureName(filename, substructure_name) return def modifySubstructureName(mol2file, name='MOL'): """ Replace the substructure name (subst_name) in a mol2 file. ARGUMENTS mol2file (string) - name of the mol2 file to modify name (string) - new substructure name NOTES This is useful becuase the OpenEye tools leave this name set to <0>. The transformation is only applied to the first molecule in the mol2 file. TODO This function is still difficult to read. It should be rewritten to be comprehensible by humans. """ # Read mol2 file. file = open(mol2file, 'r') text = file.readlines() file.close() # Find the atom records. atomsec = list() ct = 0 while text[ct].find('ATOM') == -1: ct += 1 ct += 1 atomstart = ct while text[ct].find('BOND') == -1: ct += 1 atomend = ct atomsec = text[atomstart:atomend] outtext=text[0:atomstart] repltext = atomsec[0].split()[7] # mol2 file uses space delimited, not fixed-width # Replace substructure name. for line in atomsec: # If we blindly search and replace, we'll tend to clobber stuff, as the subst_name might be "1" or something lame like that that will occur all over. # If it only occurs once, just replace it. if line.count(repltext)==1: outtext.append( line.replace(repltext, name) ) else: # Otherwise grab the string left and right of the subst_name and sandwich the new subst_name in between. This can probably be done easier in Python 2.5 with partition, but 2.4 is still used someplaces. # Loop through the line and tag locations of every non-space entry blockstart=[] ct=0 c=' ' for ct in range(len(line)): lastc = c c = line[ct] if lastc.isspace() and not c.isspace(): blockstart.append(ct) line = line[0:blockstart[7]] + line[blockstart[7]:].replace(repltext, name, 1) outtext.append(line) # Append rest of file. for line in text[atomend:]: outtext.append(line) # Write out modified mol2 file, overwriting old one. file = open(mol2file,'w') file.writelines(outtext) file.close() return def formalCharge(molecule): """ Report the net formal charge of a molecule. ARGUMENTS molecule (OEMol) - the molecule whose formal charge is to be determined RETURN VALUES formal_charge (integer) - the net formal charge of the molecule EXAMPLE net_charge = formalCharge(molecule) """ # Create a deep copy of the molecule. molecule_copy = openeye.oechem.OEMol(molecule) # Assign formal charges. openeye.oechem.OEFormalPartialCharges(molecule_copy) # Compute net formal charge. formal_charge = int(round(openeye.oechem.OENetCharge(molecule_copy))) # return formal charge return formal_charge def normalizeMolecule(molecule): """ Normalize the molecule by checking aromaticity, adding explicit hydrogens, and renaming by IUPAC name. ARGUMENTS molecule (OEMol) - the molecule to be normalized. EXAMPLES >>> # read a partial molecule and normalize it >>> molecule = readMolecule('molecule.sdf') >>> normalizeMolecule(molecule) """ # Find ring atoms and bonds openeye.oechem.OEFindRingAtomsAndBonds(molecule) # Assign aromaticity. openeye.oechem.OEAssignAromaticFlags(molecule, openeye.oechem.OEAroModelOpenEye) # Add hydrogens. openeye.oechem.OEAddExplicitHydrogens(molecule) # Set title to IUPAC name. name = openeye.oeiupac.OECreateIUPACName(molecule) molecule.SetTitle(name) return molecule def createMoleculeFromIUPAC(name, verbose=False, charge=None): """ Generate a small molecule from its IUPAC name. ARGUMENTS name (string) - IUPAC name of molecule to generate OPTIONAL ARGUMENTS verbose (boolean) - if True, subprocess output is shown (default: False) charge (int) - if specified, a form of this molecule with the desired charge state will be produced (default: None) RETURNS molecule (OEMol) - the molecule NOTES OpenEye LexiChem's OEParseIUPACName is used to generate the molecle. The molecule is normalized by adding hydrogens. Omega is used to generate a single conformation. Also note that atom names will be blank coming from this molecule. They are assigned when the molecule is written, or one can assign using OETriposAtomNames for example. EXAMPLES >>> # Generate a mol2 file for phenol. >>> molecule = createMoleculeFromIUPAC('phenol') """ # Create an OEMol molecule from IUPAC name. molecule = openeye.oechem.OEMol() # create a molecule status = openeye.oeiupac.OEParseIUPACName(molecule, name) # populate the molecule from the IUPAC name # Normalize the molecule. #normalizeMolecule(molecule) # Generate a conformation with Omega. omega = openeye.oeomega.OEOmega() omega.SetIncludeInput(False) # don't include input omega.SetMaxConfs(1) # set maximum number of conformations to 1 omega(molecule) # generate conformation if (charge != None): # Enumerate protonation states and select desired state. protonation_states = enumerateStates(molecule, enumerate="protonation", verbose=verbose) for molecule in protonation_states: if formalCharge(molecule) == charge: # Return the molecule if we've found one in the desired protonation state. return molecule if formalCharge(molecule) != charge: print "enumerateStates did not enumerate a molecule with desired formal charge." print "Options are:" for molecule in protonation_states: print "%s, formal charge %d" % (molecule.GetTitle(), formalCharge(molecule)) raise "Could not find desired formal charge." # Return the molecule. return molecule def antechamber(molecule, charge_model=False, judgetypes=None, cleanup=True, show_warnings=True, verbose=False, resname=None, netcharge=None, mm=None): """ Parameterize small molecule with GAFF using AmberTools 'antechamber' and create an OpenMM System object. ARGUMENTS molecule (OEMol) - molecule to parameterize (only the first configuration will be used if multiple are present) OPTIONAL ARGUMENTS charge_model (string) - if not False, antechamber is used to assign charges (default: False) -- if set to 'bcc', for example, AM1-BCC charges will be used judgetypes (integer) - if provided, argument passed to -j of antechamber to judge types (default: None) cleanup (boolean) - clean up temporary files (default: True) show_warnings (boolean) - show warnings during parameterization (default: True) verbose (boolean) - show all output from subprocesses (default: False) resname (string) - if set, residue name to use for parameterized molecule (default: None) netcharge (integer) -- if set, pass this net charge to calculation in antechamber (with -nc (netcharge)), otherwise assumes zero. mm (simtk.chem.openmm or equivalent) - OpenMM API to use (default: simtk.chem.openmm) RETURNS system (mm.System object) - the parameterized small molecule coordinates (simtk.unit.Quantity of natoms x 3 numpy array) - coordinates REQUIREMENTS AmberTools installation (in PATH) EXAMPLES >>> # create a molecule >>> molecule = createMoleculeFromIUPAC('phenol') >>> # parameterize it for AMBER, using antechamber to assign AM1-BCC charges >>> [system, coordinates] = antechamber(molecule, charge_model='bcc') """ # Defaults if mm is None: mm = simtk.chem.openmm # Create temporary working directory and copy ligand mol2 file there. working_directory = tempfile.mkdtemp() old_directory = os.getcwd() os.chdir(working_directory) if verbose: print "Working directory is %(working_directory)s" % vars() # Write molecule to mol2 file. tripos_mol2_filename = os.path.join(working_directory, 'tripos.mol2') writeMolecule(molecule, tripos_mol2_filename) if resname is not None: # Set substructure name (which will become residue name) if desired. modifySubstructureName(tripos_mol2_filename, resname) # Run antechamber to assign GAFF atom types. gaff_mol2_filename = os.path.join(working_directory, 'gaff.mol2') if netcharge: chargstr = '-nc %d' % netcharge else: chargestr='' command = 'antechamber -i %(tripos_mol2_filename)s -fi mol2 -o %(gaff_mol2_filename)s -fo mol2 %(chargestr)s' % vars() if judgetypes: command += ' -j %(judgetypes)d' % vars() if charge_model: formal_charge = formalCharge(molecule) command += ' -c %(charge_model)s -nc %(formal_charge)d' % vars() if verbose: print command output = commands.getoutput(command) if verbose or (output.find('Warning')>=0): print output # Generate frcmod file for additional GAFF parameters. frcmod_filename = os.path.join(working_directory, 'gaff.frcmod') print commands.getoutput('parmchk -i %(gaff_mol2_filename)s -f mol2 -o %(frcmod_filename)s' % vars()) # Create AMBER topology/coordinate files using LEaP. leapscript = """\ source leaprc.gaff mods = loadAmberParams %(frcmod_filename)s molecule = loadMol2 %(gaff_mol2_filename)s desc molecule check molecule saveAmberParm molecule amber.prmtop amber.crd quit""" % vars() leapin_filename = os.path.join(working_directory, 'leap.in') outfile = open(leapin_filename, 'w') outfile.write(leapscript) outfile.close() tleapout = commands.getoutput('tleap -f %(leapin_filename)s' % vars()) if verbose: print tleapout tleapout = tleapout.split('\n') # Shop any warnings. if show_warnings: print "Any LEaP warnings follow:" # TODO: Only print this if there are any warnings to be printed. for line in tleapout: tmp = line.upper() if tmp.find('WARNING')>-1: print line if tmp.find('ERROR')>-1: print line # Restore old directory. os.chdir(old_directory) # Create OpenMM system and coordinates from reading AMBER prmtop and crd files. import simtk.chem.openmm.extras.amber as amber prmtop_filename = os.path.join(working_directory, 'amber.prmtop') system = amber.readAmberSystem(prmtop_filename, shake=None, gbmodel=None, nonbondedCutoff=None, nonbondedMethod='no-cutoff', scee=1.2, scnb=2.0, mm=mm, verbose=verbose, EwaldErrorTolerance=None, flexibleConstraints=True) crd_filename = os.path.join(working_directory, 'amber.crd') coordinates = amber.readAmberCoordinates(crd_filename, verbose=verbose) # Clean up temporary files. os.chdir(old_directory) if cleanup: commands.getoutput('rm -r %s' % working_directory) else: print "Work done in %s..." % working_directory return [system, coordinates] #============================================================================================= # MAIN AND TESTS #============================================================================================= if __name__ == "__main__": import doctest doctest.testmod()