Merge pull request #93 from kyleabeauchamp/checkbool

Raise exception if assign charges fails

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)

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)