diff --git a/gaff2xml/openeye.py b/gaff2xml/openeye.py index bbf43ac..48fbfc6 100644 --- a/gaff2xml/openeye.py +++ b/gaff2xml/openeye.py @@ -21,7 +21,8 @@ def get_charges(molecule, max_confs=800, strictStereo=True, keep_confs=None): max_confs : int, optional, default=800 Max number of conformers to generate strictStereo : bool, optional, default=True - Adhere to strict specification of stereo isomer + If False, permits smiles strings with unspecified stereochemistry. + See https://docs.eyesopen.com/omega/usage.html keep_confs : int, optional, default=None If not None, only keep this many conformations in the final charged OEMol. Multiple conformations are still used to *determine* @@ -48,7 +49,10 @@ def get_charges(molecule, max_confs=800, strictStereo=True, keep_confs=None): charged_copy = generate_conformers(molecule, max_confs=max_confs, strictStereo=strictStereo) # Generate up to max_confs conformers - oequacpac.OEAssignPartialCharges(charged_copy, oequacpac.OECharges_AM1BCCSym) # AM1BCCSym recommended by Chris Bayly to KAB+JDC, Oct. 20 2014. + status = oequacpac.OEAssignPartialCharges(charged_copy, oequacpac.OECharges_AM1BCCSym) # AM1BCCSym recommended by Chris Bayly to KAB+JDC, Oct. 20 2014. + + if not status: + raise(RuntimeError("OEAssignPartialCharges returned error code %d" % status)) for k, conf in enumerate(charged_copy.GetConfs()): if keep_confs is not None and k > keep_confs - 1: @@ -160,7 +164,7 @@ def generate_conformers(molecule, max_confs=800, strictStereo=True, ewindow=15.0 max_confs : int, optional, default=800 Max number of conformers to generate. If None, use default OE Value. strictStereo : bool, optional, default=True - Adhere to strict specification of stereo isomer + If False, permits smiles strings with unspecified stereochemistry. Returns ------- @@ -196,7 +200,12 @@ def generate_conformers(molecule, max_confs=800, strictStereo=True, ewindow=15.0 omega.SetIncludeInput(False) # don't include input if max_confs is not None: omega.SetMaxConfs(max_confs) - omega(molcopy) # generate conformation + + status = omega(molcopy) # generate conformation + + if not status: + raise(RuntimeError("omega returned error code %d" % status)) + return molcopy @@ -346,7 +355,7 @@ def oemols_to_ffxml(molecules, base_molecule_name="lig"): return all_trajectories, ffxml -def smiles_to_antechamber(smiles_string, gaff_mol2_filename, frcmod_filename, residue_name="MOL"): +def smiles_to_antechamber(smiles_string, gaff_mol2_filename, frcmod_filename, residue_name="MOL", strictStereo=False): """Build a molecule from a smiles string and run antechamber, generating GAFF mol2 and frcmod files from a smiles string. Charges will be generated using the OpenEye QuacPac AM1-BCC implementation. @@ -366,16 +375,23 @@ def smiles_to_antechamber(smiles_string, gaff_mol2_filename, frcmod_filename, re This chokes many mol2 parsers, so we replace it with a string of your choosing. This might be useful for downstream applications if the residue names are required to be unique. + strictStereo : bool, optional, default=False + If False, permits smiles strings with unspecified stereochemistry. + See https://docs.eyesopen.com/omega/usage.html """ + oechem = import_("openeye.oechem") + if not oechem.OEChemIsLicensed(): raise(ImportError("Need License for oechem!")) + # Get the absolute path so we can find these filenames from inside a temporary directory. gaff_mol2_filename = os.path.abspath(gaff_mol2_filename) frcmod_filename = os.path.abspath(frcmod_filename) m = smiles_to_oemol(smiles_string) - m = get_charges(m, strictStereo=False, keep_confs=1) + m = get_charges(m, strictStereo=strictStereo, keep_confs=1) with enter_temp_directory(): # Avoid dumping 50 antechamber files in local directory. _unused = molecule_to_mol2(m, "./tmp.mol2", residue_name=residue_name) - tmp_gaff_mol2_filename, tmp_frcmod_filename = run_antechamber("tmp", "./tmp.mol2", charge_method=None) # USE OE AM1BCC charges! + net_charge = oechem.OENetCharge(m) + tmp_gaff_mol2_filename, tmp_frcmod_filename = run_antechamber("tmp", "./tmp.mol2", charge_method=None, net_charge=net_charge) # USE OE AM1BCC charges! shutil.copy(tmp_gaff_mol2_filename, gaff_mol2_filename) shutil.copy(tmp_frcmod_filename, frcmod_filename) diff --git a/gaff2xml/tests/test_openeye.py b/gaff2xml/tests/test_openeye.py index 604902d..45b4b01 100644 --- a/gaff2xml/tests/test_openeye.py +++ b/gaff2xml/tests/test_openeye.py @@ -9,6 +9,9 @@ import gaff2xml.openeye import pandas as pd import mdtraj as md +from mdtraj.testing import raises + +smiles_fails_with_strictStereo = "CN1CCN(CC1)CCCOc2cc3c(cc2OC)C(=[NH+]c4cc(c(cc4Cl)Cl)OC)C(=C=[N-])C=[NH+]3" try: oechem = utils.import_("openeye.oechem") @@ -183,3 +186,25 @@ def test_ffxml_simulation(): simulation.context.setPositions(model.positions) print("running") simulation.step(1) + +@skipIf(not HAVE_OE, "Cannot test openeye module without OpenEye tools.") +@raises(RuntimeError) +def test_charge_fail1(): + with utils.enter_temp_directory(): + gaff2xml.openeye.smiles_to_antechamber(smiles_fails_with_strictStereo, "test.mol2", "test.frcmod", strictStereo=True) + +@skipIf(not HAVE_OE, "Cannot test openeye module without OpenEye tools.") +@raises(RuntimeError) +def test_charge_fail2(): + m = gaff2xml.openeye.smiles_to_oemol(smiles_fails_with_strictStereo) + m = gaff2xml.openeye.get_charges(m, strictStereo=True, keep_confs=1) + +@skipIf(not HAVE_OE, "Cannot test openeye module without OpenEye tools.") +def test_charge_success1(): + with utils.enter_temp_directory(): + gaff2xml.openeye.smiles_to_antechamber(smiles_fails_with_strictStereo, "test.mol2", "test.frcmod", strictStereo=False) + +@skipIf(not HAVE_OE, "Cannot test openeye module without OpenEye tools.") +def test_charge_success2(): + m = gaff2xml.openeye.smiles_to_oemol(smiles_fails_with_strictStereo) + m = gaff2xml.openeye.get_charges(m, strictStereo=False)