Skip to content

Commit

Permalink
Finalise test and handling for multi-conformer (aka trajectory) conta…
Browse files Browse the repository at this point in the history
…ining OBMols
  • Loading branch information
lunamorrow committed Sep 10, 2024
1 parent 2ae74d3 commit 491a7f9
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 28 deletions.
29 changes: 21 additions & 8 deletions mda_openbabel_converter/OpenBabel.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,15 +46,28 @@ def __init__(self, filename: OBMol, **kwargs):
"""
Converts file to OBMol to AtomGroup
"""
n_conf = filename.NumConformers()
if n_conf > 1:
raise("Reader does not support OBMol's with more than one conformer")
n_atoms = filename.NumAtoms()
coordinates = np.array([
[(coords := atom.GetVector()).GetX(),
coords.GetY(),
coords.GetZ()] for atom in OBMolAtomIter(filename)],
dtype=np.float32)
# single position
if filename.NumConformers() == 1:
coordinates = np.array([
[(coords := atom.GetVector()).GetX(),
coords.GetY(),
coords.GetZ()] for atom in OBMolAtomIter(filename)],
dtype=np.float32)
else:
# multiple conformers, such as for a trajectory
numConf = filename.NumConformers()
coordinates = np.zeros((numConf, n_atoms, 3))
for conf_id in range(numConf):
filename.SetConformer(conf_id)
for atom in OBMolAtomIter(filename):
coordinates_inner = np.array([
[(coords := atom.GetVector()).GetX(),
coords.GetY(),
coords.GetZ()] for atom in OBMolAtomIter(filename)],
dtype=np.float32)
coordinates[conf_id] = coordinates_inner
# no coordinates present
if not np.any(coordinates):
warnings.warn("No coordinates found in the OBMol")
coordinates = np.empty((1, n_atoms, 3), dtype=np.float32)
Expand Down
44 changes: 24 additions & 20 deletions mda_openbabel_converter/tests/test_openbabel_reader.py
Original file line number Diff line number Diff line change
Expand Up @@ -67,27 +67,31 @@ def test_coordinates_complex(self, n_frames):
dtype=np.float32)
assert_equal(expected, universe.trajectory.coordinate_array[0])

# def test_multiple_conformer(self):
# obconversion = ob.OBConversion()
# obconversion.SetInFormat("xyz")
# obmol = ob.OBMol()
# notatend = obconversion.ReadFile(obmol, MET_ENKAPH_MOVIE.as_posix())
# obmol_master = ob.OBMol()
# obconversion.ReadFile(obmol_master, MET_ENKAPH_MOVIE.as_posix())
# index = 0
# while notatend:
# # temp = obmol_master.GetConformer(0)
# # obmol.CopyConformer(temp, 0)
# #obmol_master.AddConformer(temp)
def test_multiple_conformer(self):
ob_conversion = ob.OBConversion()
ob_conversion.SetInFormat("smi")
mol = ob.OBMol()
ob_conversion.ReadString(mol, "CC(=O)OC1=CC=CC=C1C(=O)O")
assert mol.AddHydrogens()

# print(obmol.GetConformer(0))
# temp = obmol.GetConformer(0)
# obmol_master.AddConformer(temp)

# obmol = ob.OBMol()
# notatend = obconversion.Read(obmol)
# print(obmol_master.NumConformers())
# assert 1 == 0
# build basic 3D coordinates
builder = ob.OBBuilder()
assert builder.Build(mol)

# generate multiconfs
confsearch = ob.OBConformerSearch()
num_conformers = 3
confsearch.Setup(mol, num_conformers)
confsearch.Search()
confsearch.GetConformers(mol)
assert mol.NumConformers() == 3

# convert to Universe
universe = mda.Universe(mol)
n_atoms = mol.NumAtoms()
expected_coordinates = np.empty((mol.NumConformers(), n_atoms, 3))
expected_shape = expected_coordinates.shape
assert_equal(universe.trajectory.coordinate_array.shape, expected_shape)

def test_no_coordinates(self):
obConversion = ob.OBConversion()
Expand Down

0 comments on commit 491a7f9

Please sign in to comment.