From 0d29212c5173033a905813935a58346c4323dcbe Mon Sep 17 00:00:00 2001 From: j-wags Date: Fri, 9 Nov 2018 13:04:43 -0800 Subject: [PATCH] Added core properties to test_to_from_rdkit. Partway through implementing/debugging stereochemistry in RDKit conversions. --- openforcefield/tests/test_molecule.py | 19 +- openforcefield/tests/test_toolkits.py | 304 +++++++++++++++++--------- openforcefield/topology/molecule.py | 6 +- openforcefield/utils/toolkits.py | 28 ++- 4 files changed, 234 insertions(+), 123 deletions(-) diff --git a/openforcefield/tests/test_molecule.py b/openforcefield/tests/test_molecule.py index b90538b32..daa24ecab 100644 --- a/openforcefield/tests/test_molecule.py +++ b/openforcefield/tests/test_molecule.py @@ -584,6 +584,7 @@ def test_name(self): molecule.name = name assert molecule.name == name + # TODO: This should be a toolkit test #@pytest.mark.skipif(OPENEYE_UNAVAILABLE, reason=_OPENEYE_UNAVAILABLE_MESSAGE) @OpenEyeToolkitWrapper.requires_toolkit() def test_iupac_roundtrip(self): @@ -600,6 +601,7 @@ def test_topology_roundtrip(self): molecule_copy = Molecule.from_topology(topology) assert molecule == molecule_copy + # TODO: This should be a toolkit test def test_file_roundtrip(self): """Test to/from file""" import os @@ -626,15 +628,16 @@ def test_file_roundtrip(self): # NOTE: We can't read pdb files and expect chemical information to be preserved os.unlink(iofile.name) + # Really big round-trip tests have now been added to test_toolkits and aren't necessary here #@pytest.mark.skipif(RDKIT_UNAVAILABLE, reason=_RDKIT_UNAVAILABLE_MESSAGE) - @RDKitToolkitWrapper.requires_toolkit() - def test_rdkit_roundtrip(self): - for molecule in self.molecules: - rdmol = molecule.to_rdkit() - molecule2 = Molecule.from_rdkit(rdmol) - assert_molecule_is_equal(molecule, molecule2, "Molecule.to_rdkit()/from_rdkit() round trip failed") - molecule3 = Molecule(rdmol) - assert_molecule_is_equal(molecule, molecule3, "Molecule(rdmol) constructor failed") + ##@RDKitToolkitWrapper.requires_toolkit() + #def test_rdkit_roundtrip(self): + # for molecule in self.molecules: + # rdmol = molecule.to_rdkit() + # molecule2 = Molecule.from_rdkit(rdmol) + 3 assert_molecule_is_equal(molecule, molecule2, "Molecule.to_rdkit()/from_rdkit() round trip failed") + # molecule3 = Molecule(rdmol) + # assert_molecule_is_equal(molecule, molecule3, "Molecule(rdmol) constructor failed") #@pytest.mark.skipif(OPENEYE_UNAVAILABLE, reason=_OPENEYE_UNAVAILABLE_MESSAGE) @OpenEyeToolkitWrapper.requires_toolkit() diff --git a/openforcefield/tests/test_toolkits.py b/openforcefield/tests/test_toolkits.py index c4eb0ec17..f60e82542 100644 --- a/openforcefield/tests/test_toolkits.py +++ b/openforcefield/tests/test_toolkits.py @@ -91,73 +91,79 @@ def test_smiles_charged(self): smiles2 = molecule.to_smiles(toolkit_registry=toolkit_wrapper) assert smiles == smiles2 - - @pytest.mark.skipif( not OpenEyeToolkitWrapper.toolkit_is_available(), reason='OpenEye Toolkit not available') - def test_to_from_openeye_core_props_filled(self): - """Test OpenEyeToolkitWrapper to_openeye() and from_openeye()""" - toolkit_wrapper = OpenEyeToolkitWrapper() - - # Replacing with a simple molecule with stereochemistry - input_smiles = 'C\C(F)=C(/F)C[C@@](C)(Cl)Br' - expected_output_smiles = '[H]C([H])([H])C(=C(C([H])([H])C(C([H])([H])[H])(Cl)Br)F)F' - molecule = Molecule.from_smiles(input_smiles, toolkit_registry=toolkit_wrapper) - assert molecule.to_smiles(toolkit_registry=toolkit_wrapper) == expected_output_smiles - - # Populate core molecule property fields - molecule.name = 'Alice' - partial_charges = unit.Quantity(np.array([-.9, -.8, -.7, -.6, - -.5, -.4, -.3, -.2, - -.1, 0., .1, .2, - .3, .4, .5, .6, - .7, .8]), unit.elementary_charge) - molecule.set_partial_charges(partial_charges) - coords = unit.Quantity(np.array([['0.0', '1.0', '2.0'], ['3.0', '4.0', '5.0'], ['6.0', '7.0', '8.0'], - ['9.0', '10.0', '11.0'] , ['12.0', '13.0', '14.0'], ['15.0', '16.0', '17.0'], - ['18.0', '19.0', '20.0'], ['21.0', '22.0', '23.0'], ['24.0', '25.0', '26.0'], - ['27.0', '28.0', '29.0'], ['30.0', '31.0', '32.0'], ['33.0', '34.0', '35.0'], - ['36.0', '37.0', '38.0'], ['39.0', '40.0', '41.0'], ['42.0', '43.0', '44.0'], - ['45.0', '46.0', '47.0'], ['48.0', '49.0', '50.0'], ['51.0', '52.0', '53.0']]), - unit.angstrom) - molecule.add_conformer(coords) - # Populate core atom property fields - molecule.atoms[2].name = 'Bob' - # Ensure one atom has its stereochemistry specified - central_carbon_stereo_specified = False - for atom in molecule.atoms: - if (atom.atomic_number == 6) and atom.stereochemistry == "S": - central_carbon_stereo_specified = True - assert central_carbon_stereo_specified - - # Populate bond core property fields - fractional_bond_orders = [float(val) for val in range(18)] - for fbo, bond in zip(fractional_bond_orders, molecule.bonds): - bond.fractional_bond_order = fbo - - # Do a first conversion to/from oemol - oemol = molecule.to_openeye() - molecule2 = Molecule.from_openeye(oemol) - - # Test that properties survived first conversion - #assert molecule.to_dict() == molecule2.to_dict() - assert molecule.name == molecule2.name - # NOTE: This expects the same indexing scheme in the original and new molecule - - central_carbon_stereo_specified = False - for atom in molecule2.atoms: - if (atom.atomic_number == 6) and atom.stereochemistry == "S": - central_carbon_stereo_specified = True - assert central_carbon_stereo_specified - for atom1, atom2 in zip(molecule.atoms, molecule2.atoms): - assert atom1.to_dict() == atom2.to_dict() - for bond1, bond2 in zip(molecule.bonds, molecule2.bonds): - assert bond1.to_dict() == bond2.to_dict() - assert (molecule._conformers[0] == molecule2._conformers[0]).all() - for pc1, pc2 in zip(molecule._partial_charges, molecule2._partial_charges): - pc1_ul = pc1 / unit.elementary_charge - pc2_ul = pc2 / unit.elementary_charge - assert_almost_equal(pc1_ul, pc2_ul, decimal=6) - assert molecule2.to_smiles(toolkit_registry=toolkit_wrapper) == expected_output_smiles - # TODO: This should be its own test + # TODO: Make a close copy of this test, but with core props unset, and make sure they are equal to the default value both before and after toolkit conversion + + @pytest.mark.skipif(not OpenEyeToolkitWrapper.toolkit_is_available(), reason='OpenEye Toolkit not available') + def test_to_from_openeye_core_props_filled(self): + """Test OpenEyeToolkitWrapper to_openeye() and from_openeye()""" + toolkit_wrapper = OpenEyeToolkitWrapper() + + # Replacing with a simple molecule with stereochemistry + input_smiles = 'C\C(F)=C(/F)C[C@@](C)(Cl)Br' + expected_output_smiles = '[H]C([H])([H])C(=C(C([H])([H])C(C([H])([H])[H])(Cl)Br)F)F' + molecule = Molecule.from_smiles(input_smiles, toolkit_registry=toolkit_wrapper) + assert molecule.to_smiles(toolkit_registry=toolkit_wrapper) == expected_output_smiles + + # Populate core molecule property fields + molecule.name = 'Alice' + partial_charges = unit.Quantity(np.array([-.9, -.8, -.7, -.6, + -.5, -.4, -.3, -.2, + -.1, 0., .1, .2, + .3, .4, .5, .6, + .7, .8]), unit.elementary_charge) + molecule.set_partial_charges(partial_charges) + coords = unit.Quantity(np.array([['0.0', '1.0', '2.0'], ['3.0', '4.0', '5.0'], ['6.0', '7.0', '8.0'], + ['9.0', '10.0', '11.0'], ['12.0', '13.0', '14.0'], + ['15.0', '16.0', '17.0'], + ['18.0', '19.0', '20.0'], ['21.0', '22.0', '23.0'], + ['24.0', '25.0', '26.0'], + ['27.0', '28.0', '29.0'], ['30.0', '31.0', '32.0'], + ['33.0', '34.0', '35.0'], + ['36.0', '37.0', '38.0'], ['39.0', '40.0', '41.0'], + ['42.0', '43.0', '44.0'], + ['45.0', '46.0', '47.0'], ['48.0', '49.0', '50.0'], + ['51.0', '52.0', '53.0']]), + unit.angstrom) + molecule.add_conformer(coords) + # Populate core atom property fields + molecule.atoms[2].name = 'Bob' + # Ensure one atom has its stereochemistry specified + central_carbon_stereo_specified = False + for atom in molecule.atoms: + if (atom.atomic_number == 6) and atom.stereochemistry == "S": + central_carbon_stereo_specified = True + assert central_carbon_stereo_specified + + # Populate bond core property fields + fractional_bond_orders = [float(val) for val in range(18)] + for fbo, bond in zip(fractional_bond_orders, molecule.bonds): + bond.fractional_bond_order = fbo + + # Do a first conversion to/from oemol + oemol = molecule.to_openeye() + molecule2 = Molecule.from_openeye(oemol) + + # Test that properties survived first conversion + # assert molecule.to_dict() == molecule2.to_dict() + assert molecule.name == molecule2.name + # NOTE: This expects the same indexing scheme in the original and new molecule + + central_carbon_stereo_specified = False + for atom in molecule2.atoms: + if (atom.atomic_number == 6) and atom.stereochemistry == "S": + central_carbon_stereo_specified = True + assert central_carbon_stereo_specified + for atom1, atom2 in zip(molecule.atoms, molecule2.atoms): + assert atom1.to_dict() == atom2.to_dict() + for bond1, bond2 in zip(molecule.bonds, molecule2.bonds): + assert bond1.to_dict() == bond2.to_dict() + assert (molecule._conformers[0] == molecule2._conformers[0]).all() + for pc1, pc2 in zip(molecule._partial_charges, molecule2._partial_charges): + pc1_ul = pc1 / unit.elementary_charge + pc2_ul = pc2 / unit.elementary_charge + assert_almost_equal(pc1_ul, pc2_ul, decimal=6) + assert molecule2.to_smiles(toolkit_registry=toolkit_wrapper) == expected_output_smiles + # TODO: This should be its own test #charges = unit.Quantity(np.array([ -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35], dtype=np.float), # unit.elementary_charge) @@ -234,6 +240,9 @@ def test_generate_conformers(self): smiles = '[H]C([H])([H])C([H])([H])[H]' molecule = toolkit_wrapper.from_smiles(smiles) molecule.generate_conformers() + assert molecule.n_conformers != 0 + assert not(molecule.conformers[0] == (0.*unit.angstrom)).all() + # TODO: Make this test more robust @pytest.mark.skipif( not OpenEyeToolkitWrapper.toolkit_is_available(), reason='OpenEye Toolkit not available') @@ -303,8 +312,7 @@ def test_compute_wiberg_bond_orders(self): print([bond.fractional_bond_order for bond in molecule.bonds]) # TODO: Add test for equivalent Wiberg orders for equivalent bonds - - @pytest.mark.skipif( not OpenEyeToolkitWrapper.toolkit_is_available(), reason='OpenEye Toolkit not available') + @pytest.mark.skipif(not OpenEyeToolkitWrapper.toolkit_is_available(), reason='OpenEye Toolkit not available') def test_compute_wiberg_bond_orders_charged(self): """Test OpenEyeToolkitWrapper compute_wiberg_bond_orders() on a molecule with net charge +1""" @@ -312,18 +320,36 @@ def test_compute_wiberg_bond_orders_charged(self): smiles = '[H]C([H])([H])[N+]([H])([H])[H]' molecule = toolkit_wrapper.from_smiles(smiles) molecule.generate_conformers(toolkit_registry=toolkit_wrapper) - for charge_model in ['am1','pm3']: + for charge_model in ['am1', 'pm3']: molecule.compute_wiberg_bond_orders(toolkit_registry=toolkit_wrapper, charge_model=charge_model) # TODO: Add test for equivalent Wiberg orders for equivalent bonds + @pytest.mark.skipif(not OpenEyeToolkitWrapper.toolkit_is_available(), + reason='OpenEye Toolkit not available') + def test_compute_wiberg_bond_orders_double_bond(self): + """Test OpenEyeToolkitWrapper compute_wiberg_bond_orders() on a molecule with net charge +1""" + + toolkit_wrapper = OpenEyeToolkitWrapper() + smiles = 'C\C(F)=C(/F)C[C@@](C)(Cl)Br' + molecule = toolkit_wrapper.from_smiles(smiles) + molecule.generate_conformers(toolkit_registry=toolkit_wrapper) + for charge_model in ['am1', 'pm3']: + molecule.compute_wiberg_bond_orders(toolkit_registry=toolkit_wrapper, charge_model=charge_model) + # TODO: Add test for equivalent Wiberg orders for equivalent bonds + + double_bond_has_wbo_near_2 = False + for bond in molecule.bonds: + if bond.bond_order == 2: + if 1.75 < bond.fractional_bond_order < 2.25: + double_bond_has_wbo_near_2 = True + assert double_bond_has_wbo_near_2 + # TODO: Check partial charge invariants (total charge, charge equivalence) - # TODO: Add test for partial charge calculation with formally charged atoms/molecules - - # TODO: Add test for higher bonds orders + # TODO: Add test for aromaticity # TODO: Add test and molecule functionality for isotopes @@ -370,39 +396,111 @@ def test_smiles_charged(self): smiles2 = molecule.to_smiles(toolkit_registry=toolkit_wrapper) assert smiles == smiles2 - @pytest.mark.skipif( not RDKitToolkitWrapper.toolkit_is_available(), reason='RDKit Toolkit not available') - def test_to_from_rdkit(self): - """Test RDKitToolkitWrapper to_rdkit() and from_rdkit()""" + + + # TODO: Make a close copy of this test, but with core props unset, and make sure they are equal to the default value both before and after toolkit conversion + + @pytest.mark.skipif( not RDKitToolkitWrapper.toolkit_is_available(), reason='OpenEye Toolkit not available') + def test_to_from_rdkit_core_props_filled(self): + """Test OpenEyeToolkitWrapper to_openeye() and from_openeye()""" toolkit_wrapper = RDKitToolkitWrapper() - # This differs from OE's expected output due to different canonicalization schemes - smiles = '[H][C]([H])([H])[C]([H])([H])[H]' - molecule = Molecule.from_smiles(smiles, - toolkit_registry=toolkit_wrapper) + + # Replacing with a simple molecule with stereochemistry + input_smiles = 'C\C(F)=C(/F)C[C@@](C)(Cl)Br' + expected_output_smiles = '[H]C([H])([H])C(=C(C([H])([H])C(C([H])([H])[H])(Cl)Br)F)F' + molecule = Molecule.from_smiles(input_smiles, toolkit_registry=toolkit_wrapper) + assert molecule.to_smiles(toolkit_registry=toolkit_wrapper) == expected_output_smiles + + # Populate core molecule property fields + molecule.name = 'Alice' + partial_charges = unit.Quantity(np.array([-.9, -.8, -.7, -.6, + -.5, -.4, -.3, -.2, + -.1, 0., .1, .2, + .3, .4, .5, .6, + .7, .8]), unit.elementary_charge) + molecule.set_partial_charges(partial_charges) + coords = unit.Quantity(np.array([['0.0', '1.0', '2.0'], ['3.0', '4.0', '5.0'], ['6.0', '7.0', '8.0'], + ['9.0', '10.0', '11.0'] , ['12.0', '13.0', '14.0'], ['15.0', '16.0', '17.0'], + ['18.0', '19.0', '20.0'], ['21.0', '22.0', '23.0'], ['24.0', '25.0', '26.0'], + ['27.0', '28.0', '29.0'], ['30.0', '31.0', '32.0'], ['33.0', '34.0', '35.0'], + ['36.0', '37.0', '38.0'], ['39.0', '40.0', '41.0'], ['42.0', '43.0', '44.0'], + ['45.0', '46.0', '47.0'], ['48.0', '49.0', '50.0'], ['51.0', '52.0', '53.0']]), + unit.angstrom) + molecule.add_conformer(coords) + # Populate core atom property fields + molecule.atoms[2].name = 'Bob' + # Ensure one atom has its stereochemistry specified + central_carbon_stereo_specified = False + for atom in molecule.atoms: + if (atom.atomic_number == 6) and atom.stereochemistry == "S": + central_carbon_stereo_specified = True + assert central_carbon_stereo_specified + + # Populate bond core property fields + fractional_bond_orders = [float(val) for val in range(18)] + for fbo, bond in zip(fractional_bond_orders, molecule.bonds): + bond.fractional_bond_order = fbo + + # Do a first conversion to/from oemol rdmol = molecule.to_rdkit() molecule2 = Molecule.from_rdkit(rdmol) - smiles2 = molecule2.to_smiles(toolkit_registry=toolkit_wrapper) - new_conf1 = unit.Quantity(np.array([[1,2,3],[4,5,6],[7,8,9], - [10,11,12],[13,14,15],[16,17,18], - [19,20,21],[22,23,24],[25,26,27]], - dtype=np.float), - unit.angstrom) - new_conf2 = unit.Quantity(np.array([[101,102,103],[104,105,106],[107,108,109], - [110,111,112],[113,114,115],[116,117,118], - [119,120,121],[122,123,124],[125,126,127]], - dtype=np.float), - unit.angstrom) - molecule2.add_conformer(new_conf1) - molecule2.add_conformer(new_conf2) - - - smiles2 = molecule2.to_smiles(toolkit_registry=toolkit_wrapper) - assert smiles == smiles2 - oemol2 = Molecule.to_openeye(molecule2) - assert oemol2.NumConfs() == 2 - molecule3 = Molecule.from_openeye(oemol2) - assert len(molecule3._conformers) == 2 - assert (molecule2._conformers[0] == molecule3._conformers[0]).all() - assert (molecule2._conformers[1] == molecule3._conformers[1]).all() + + # Test that properties survived first conversion + #assert molecule.to_dict() == molecule2.to_dict() + assert molecule.name == molecule2.name + # NOTE: This expects the same indexing scheme in the original and new molecule + + central_carbon_stereo_specified = False + for atom in molecule2.atoms: + if (atom.atomic_number == 6) and atom.stereochemistry == "S": + central_carbon_stereo_specified = True + assert central_carbon_stereo_specified + for atom1, atom2 in zip(molecule.atoms, molecule2.atoms): + assert atom1.to_dict() == atom2.to_dict() + for bond1, bond2 in zip(molecule.bonds, molecule2.bonds): + assert bond1.to_dict() == bond2.to_dict() + assert (molecule._conformers[0] == molecule2._conformers[0]).all() + for pc1, pc2 in zip(molecule._partial_charges, molecule2._partial_charges): + pc1_ul = pc1 / unit.elementary_charge + pc2_ul = pc2 / unit.elementary_charge + assert_almost_equal(pc1_ul, pc2_ul, decimal=6) + assert molecule2.to_smiles(toolkit_registry=toolkit_wrapper) == expected_output_smiles + # TODO: This should be its own test + + + #@pytest.mark.skipif( not RDKitToolkitWrapper.toolkit_is_available(), reason='RDKit Toolkit not available') + #def test_to_from_rdkit(self): + # """Test RDKitToolkitWrapper to_rdkit() and from_rdkit()""" + # toolkit_wrapper = RDKitToolkitWrapper() + # # This differs from OE's expected output due to different canonicalization schemes + # smiles = '[H][C]([H])([H])[C]([H])([H])[H]' + # molecule = Molecule.from_smiles(smiles, + # toolkit_registry=toolkit_wrapper) + # rdmol = molecule.to_rdkit() + # molecule2 = Molecule.from_rdkit(rdmol) + # smiles2 = molecule2.to_smiles(toolkit_registry=toolkit_wrapper) + # new_conf1 = unit.Quantity(np.array([[1,2,3],[4,5,6],[7,8,9], + # [10,11,12],[13,14,15],[16,17,18], + # [19,20,21],[22,23,24]], + # dtype=np.float), + # unit.angstrom) + # new_conf2 = unit.Quantity(np.array([[101,102,103],[104,105,106],[107,108,109], + # [110,111,112],[113,114,115],[116,117,118], + # [119,120,121],[122,123,124]], + # dtype=np.float), + # unit.angstrom) + # molecule2.add_conformer(new_conf1) + # molecule2.add_conformer(new_conf2) + # + # + # smiles2 = molecule2.to_smiles(toolkit_registry=toolkit_wrapper) + # assert smiles == smiles2 + # oemol2 = Molecule.to_openeye(molecule2) + # assert oemol2.NumConfs() == 2 + # molecule3 = Molecule.from_openeye(oemol2) + # assert len(molecule3._conformers) == 2 + # assert (molecule2._conformers[0] == molecule3._conformers[0]).all() + # assert (molecule2._conformers[1] == molecule3._conformers[1]).all() @pytest.mark.skipif( not RDKitToolkitWrapper.toolkit_is_available(), reason='RDKit Toolkit not available') diff --git a/openforcefield/topology/molecule.py b/openforcefield/topology/molecule.py index e441dc7fa..f7afd2adc 100644 --- a/openforcefield/topology/molecule.py +++ b/openforcefield/topology/molecule.py @@ -372,9 +372,11 @@ def bonds(self): # yield bond @property - def bonded_to(self): + #def bonded_to(self): + def bonded_atoms(self): + """ - The list of ``Atom`` objects this atom is involved in + The list of ``Atom`` objects this atom is involved in bonds with """ for bond in self._bonds: diff --git a/openforcefield/utils/toolkits.py b/openforcefield/utils/toolkits.py index 145f0219a..19c9953e9 100644 --- a/openforcefield/utils/toolkits.py +++ b/openforcefield/utils/toolkits.py @@ -1390,10 +1390,11 @@ def from_rdkit(self, rdmol): #Chem.SanitizeMol(rdmol, Chem.SANITIZE_ALL^Chem.SANITIZE_SETAROMATICITY) #Chem.SetAromaticity(rdmol, Chem.AromaticityModel.AROMATICITY_MDL) # If RDMol has a title save it - if rdmol.HasProp("_Name"): - mol.name == rdmol.GetProp("_Name") + # TODO: This was previously "_Name" -- Is there anything special about that tag? + if rdmol.HasProp("name"): + mol.name == rdmol.GetProp("name") else: - mol.name = None + mol.name = "" # Store all properties # TODO: Should Title or _Name be a special property? @@ -1571,18 +1572,25 @@ def to_rdkit(molecule, aromaticity_model=DEFAULT_AROMATICITY_MODEL): rdbond.SetBondType(_bondtypes[bond.bond_order]) rdbond.SetIsAromatic(False) + + ## Debugging: We need to somehow clean up the molecule here to make it not crash below. Otherwise it will complain that we didn't call + #status = Chem.SanitizeMol(rdmol) + #if status == False: + # raise Exception('Unable to sanitize molecule') + #rdmol.UpdatePropertyCache() + #rdmol = rdmol.AddHs(rdmol) # Assign bond stereochemistry for bond in molecule.bonds: if bond.stereochemistry: # Determine neighbors # TODO: This API needs to be created - n1 = [n.index for n in bond.atom1.bonded_atoms if n != bond.atom2][0] - n2 = [n.index for n in bond.atom2.bonded_atoms if n != bond.atom1][0] + n1 = [n.molecule_atom_index for n in bond.atom1.bonded_atoms if n != bond.atom2][0] + n2 = [n.molecule_atom_index for n in bond.atom2.bonded_atoms if n != bond.atom1][0] # Get rdmol bonds bond_atom1_index = molecule.atoms.index(bond.atom1) bond_atom2_index = molecule.atoms.index(bond.atom2) bond1 = rdmol.GetBondBetweenAtoms(map_atoms[n1], map_atoms[bond.atom1_index]) - bond2 = rdmol.GetBondBetweenAtoms(map_atoms[bond_atom1_index], map_atoms[bon_.atom2_index]) + bond2 = rdmol.GetBondBetweenAtoms(map_atoms[bond_atom1_index], map_atoms[bond.atom2_index]) bond3 = rdmol.GetBondBetweenAtoms(map_atoms[bond_atom2_index], map_atoms[n2]) # Set arbitrary stereochemistry # Since this is relative, the first bond always goes up @@ -1814,17 +1822,17 @@ def compute_partial_charges(self, molecule, charge_model=None): # Write out molecule in SDF format ## TODO: How should we handle multiple conformers? self._rdkit_toolkit_wrapper.to_file(molecule, 'molecule.sdf', outfile_format='sdf') - os.system('ls') - os.system('cat molecule.sdf') + #os.system('ls') + #os.system('cat molecule.sdf') # Compute desired charges # TODO: Add error handling if antechamber chokes # TODO: Add something cleaner than os.system os.system("antechamber -i molecule.sdf -fi sdf -o charged.mol2 -fo mol2 -pf yes -c {} -nc {}".format(charge_model, net_charge)) - os.system('cat charged.mol2') + #os.system('cat charged.mol2') # Write out just charges os.system("antechamber -i charged.mol2 -fi mol2 -o charges2.mol2 -fo mol2 -c wc -cf charges.txt -pf yes") - os.system('cat charges.txt') + #os.system('cat charges.txt') # Check to ensure charges were actually produced if not os.path.exists('charges.txt'): # TODO: copy files into local directory to aid debugging?