Indexing Molecules¶
Molecules are collections of atoms that (at least notionally) should all be bonded to each other and represent the concept of a molecule in a chemical system.
A molecule is a molecular container for atoms, residues, chains and
segments, and is implemented via the sire.mol.Molecule class.
At a minimum, a sire.mol.Molecule contains at least one
sire.mol.Atom. There
is no requirement for this atom to be placed into a residue, chain
or segment, and so it is possible for a molecule to have zero
residues, chains or segments.
There are several classes in sire that can act as molecular containers
for molecule objects. You have already encountered one, which is the
System class, into which the molecules were first
loaded.
For example, load up the aladip system;
>>> import sire as sr
>>> mols = sr.load(sr.expand(sr.tutorial_url, ["ala.top", "ala.crd"]))
>>> print(mols)
System( name=ACE num_molecules=631 num_residues=633 num_atoms=1912 )
This has loaded 631 molecules. You can access molecules by their index in the system by simple index;
>>> print(mols[0])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
or by their name
>>> print(mols["ACE"])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
You can also index via a range
>>> print(mols[0:10])
SelectorMol( size=10
0: Molecule( ACE:2 num_atoms=22 num_residues=3 )
1: Molecule( WAT:3 num_atoms=3 num_residues=1 )
2: Molecule( WAT:4 num_atoms=3 num_residues=1 )
3: Molecule( WAT:5 num_atoms=3 num_residues=1 )
4: Molecule( WAT:6 num_atoms=3 num_residues=1 )
5: Molecule( WAT:7 num_atoms=3 num_residues=1 )
6: Molecule( WAT:8 num_atoms=3 num_residues=1 )
7: Molecule( WAT:9 num_atoms=3 num_residues=1 )
8: Molecule( WAT:10 num_atoms=3 num_residues=1 )
9: Molecule( WAT:11 num_atoms=3 num_residues=1 )
)
or via a list of indicies
>>> print(mols[ [0, 2, 4, 6] ])
SelectorMol( size=4
0: Molecule( ACE:2 num_atoms=22 num_residues=3 )
1: Molecule( WAT:4 num_atoms=3 num_residues=1 )
2: Molecule( WAT:6 num_atoms=3 num_residues=1 )
3: Molecule( WAT:8 num_atoms=3 num_residues=1 )
)
If multiple molecules are returned, then the result is held in a
SelectorMol class. This is also a molecular container,
and can be used in the same way as all other molecular containers. For example,
>>> print(mols["WAT"][0])
Molecule( WAT:3 num_atoms=3 num_residues=1 )
gives the first molecule called WAT.
You can achieve the same result using the molecule()
and molecules() functions.
>>> print(mols.molecule(0))
Molecule( ACE:2 num_atoms=22 num_residues=3 )
>>> print(mols.molecules(range(0, 10)))
SelectorMol( size=10
0: Molecule( ACE:2 num_atoms=22 num_residues=3 )
1: Molecule( WAT:3 num_atoms=3 num_residues=1 )
2: Molecule( WAT:4 num_atoms=3 num_residues=1 )
3: Molecule( WAT:5 num_atoms=3 num_residues=1 )
4: Molecule( WAT:6 num_atoms=3 num_residues=1 )
5: Molecule( WAT:7 num_atoms=3 num_residues=1 )
6: Molecule( WAT:8 num_atoms=3 num_residues=1 )
7: Molecule( WAT:9 num_atoms=3 num_residues=1 )
8: Molecule( WAT:10 num_atoms=3 num_residues=1 )
9: Molecule( WAT:11 num_atoms=3 num_residues=1 )
)
>>> print(mols.molecules("WAT").molecule(0))
Molecule( WAT:3 num_atoms=3 num_residues=1 )
Search for molecules¶
You can search for a molecule in a molecule container via their
name (molname) or number (molnum).
>>> print(mols["molname ACE"])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
>>> print(mols["molnum 2"])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
>>> print(mols["molidx 0"])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
>>> print(mols["molname WAT"])
SelectorMol( size=630
0: Molecule( WAT:3 num_atoms=3 num_residues=1 )
1: Molecule( WAT:4 num_atoms=3 num_residues=1 )
2: Molecule( WAT:5 num_atoms=3 num_residues=1 )
3: Molecule( WAT:6 num_atoms=3 num_residues=1 )
4: Molecule( WAT:7 num_atoms=3 num_residues=1 )
...
625: Molecule( WAT:628 num_atoms=3 num_residues=1 )
626: Molecule( WAT:629 num_atoms=3 num_residues=1 )
627: Molecule( WAT:630 num_atoms=3 num_residues=1 )
628: Molecule( WAT:631 num_atoms=3 num_residues=1 )
629: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
Note
Note how the name and number of the molecule are printed out,
e.g. ACE:2 is for the molecule called ACE with number 2.
You can also search by the index in the molecular container (molidx),
e.g.
>>> print(mols["molidx > 50"])
SelectorMol( size=580
0: Molecule( WAT:53 num_atoms=3 num_residues=1 )
1: Molecule( WAT:54 num_atoms=3 num_residues=1 )
2: Molecule( WAT:55 num_atoms=3 num_residues=1 )
3: Molecule( WAT:56 num_atoms=3 num_residues=1 )
4: Molecule( WAT:57 num_atoms=3 num_residues=1 )
...
575: Molecule( WAT:628 num_atoms=3 num_residues=1 )
576: Molecule( WAT:629 num_atoms=3 num_residues=1 )
577: Molecule( WAT:630 num_atoms=3 num_residues=1 )
578: Molecule( WAT:631 num_atoms=3 num_residues=1 )
579: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
Note
The molidx is the index of the molecule in its parent container.
This will depend on which container the search is performed with.
This is different to atomidx, residx etc, which are unique
identifying indicies of the atom, residue etc in their parent molecule.
You can combine the search with other identifiers, e.g.
>>> print(mols["molname WAT and molnum 3"])
Molecule( WAT:3 num_atoms=3 num_residues=1 )
and can search for multiple names or numbers
>>> print(mols["molname ACE, WAT"])
SelectorMol( size=631
0: Molecule( ACE:2 num_atoms=22 num_residues=3 )
1: Molecule( WAT:3 num_atoms=3 num_residues=1 )
2: Molecule( WAT:4 num_atoms=3 num_residues=1 )
3: Molecule( WAT:5 num_atoms=3 num_residues=1 )
4: Molecule( WAT:6 num_atoms=3 num_residues=1 )
...
626: Molecule( WAT:628 num_atoms=3 num_residues=1 )
627: Molecule( WAT:629 num_atoms=3 num_residues=1 )
628: Molecule( WAT:630 num_atoms=3 num_residues=1 )
629: Molecule( WAT:631 num_atoms=3 num_residues=1 )
630: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
>>> print(mols["molnum 5:11, 30, 40"])
SelectorMol( size=8
0: Molecule( WAT:5 num_atoms=3 num_residues=1 )
1: Molecule( WAT:6 num_atoms=3 num_residues=1 )
2: Molecule( WAT:7 num_atoms=3 num_residues=1 )
3: Molecule( WAT:8 num_atoms=3 num_residues=1 )
4: Molecule( WAT:9 num_atoms=3 num_residues=1 )
5: Molecule( WAT:10 num_atoms=3 num_residues=1 )
6: Molecule( WAT:30 num_atoms=3 num_residues=1 )
7: Molecule( WAT:40 num_atoms=3 num_residues=1 )
)
Wildcard (glob) searching is also supported for molecule names.
>>> print(mols["molname /A*/"])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
>>> print(mols["molname /wat/i"])
SelectorMol( size=630
0: Molecule( WAT:3 num_atoms=3 num_residues=1 )
1: Molecule( WAT:4 num_atoms=3 num_residues=1 )
2: Molecule( WAT:5 num_atoms=3 num_residues=1 )
3: Molecule( WAT:6 num_atoms=3 num_residues=1 )
4: Molecule( WAT:7 num_atoms=3 num_residues=1 )
...
625: Molecule( WAT:628 num_atoms=3 num_residues=1 )
626: Molecule( WAT:629 num_atoms=3 num_residues=1 )
627: Molecule( WAT:630 num_atoms=3 num_residues=1 )
628: Molecule( WAT:631 num_atoms=3 num_residues=1 )
629: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
Finding atoms, residues, chains and segments in a molecule¶
Because both System and SelectorMol
are molecular containers, they both have their own
atom(), atoms(),
residue(), residues(),
etc functions for accessing atoms, residues, chains or segments across
multiple molecules.
For example, you can get all of the atoms that are oxygens using
>>> print(mols.atoms("element O"))
SireMol::SelectorM<SireMol::Atom>( size=632
0: MolNum(2) Atom( O:6 [ 19.19, 5.44, 14.76] )
1: MolNum(2) Atom( O:16 [ 14.94, 3.17, 15.88] )
2: MolNum(3) Atom( O:23 [ 25.64, 8.50, 22.42] )
3: MolNum(4) Atom( O:26 [ 22.83, 8.93, 4.14] )
4: MolNum(5) Atom( O:29 [ 9.96, 24.84, 23.52] )
...
627: MolNum(628) Atom( O:1898 [ 22.40, 11.84, 10.08] )
628: MolNum(629) Atom( O:1901 [ 0.63, 8.69, 19.94] )
629: MolNum(630) Atom( O:1904 [ 18.69, 22.12, 9.35] )
630: MolNum(631) Atom( O:1907 [ 21.65, 7.88, 9.79] )
631: MolNum(632) Atom( O:1910 [ 9.25, 13.73, 2.29] )
)
or, even easier,
>>> print(mols["element O"])
SireMol::SelectorM<SireMol::Atom>( size=632
0: MolNum(2) Atom( O:6 [ 19.19, 5.44, 14.76] )
1: MolNum(2) Atom( O:16 [ 14.94, 3.17, 15.88] )
2: MolNum(3) Atom( O:23 [ 25.64, 8.50, 22.42] )
3: MolNum(4) Atom( O:26 [ 22.83, 8.93, 4.14] )
4: MolNum(5) Atom( O:29 [ 9.96, 24.84, 23.52] )
...
627: MolNum(628) Atom( O:1898 [ 22.40, 11.84, 10.08] )
628: MolNum(629) Atom( O:1901 [ 0.63, 8.69, 19.94] )
629: MolNum(630) Atom( O:1904 [ 18.69, 22.12, 9.35] )
630: MolNum(631) Atom( O:1907 [ 21.65, 7.88, 9.79] )
631: MolNum(632) Atom( O:1910 [ 9.25, 13.73, 2.29] )
)
The result is a SelectorM_Atom_. This is a multi-molecule
version of Selector_Atom_, which behaves in an identical
way.
This works for residues too!
>>> print(mols.residues("WAT"))
SireMol::SelectorM<SireMol::Residue>( size=630
0: MolNum(3) Residue( WAT:4 num_atoms=3 )
1: MolNum(4) Residue( WAT:5 num_atoms=3 )
2: MolNum(5) Residue( WAT:6 num_atoms=3 )
3: MolNum(6) Residue( WAT:7 num_atoms=3 )
4: MolNum(7) Residue( WAT:8 num_atoms=3 )
...
625: MolNum(628) Residue( WAT:629 num_atoms=3 )
626: MolNum(629) Residue( WAT:630 num_atoms=3 )
627: MolNum(630) Residue( WAT:631 num_atoms=3 )
628: MolNum(631) Residue( WAT:632 num_atoms=3 )
629: MolNum(632) Residue( WAT:633 num_atoms=3 )
)
>>> print(mols["resname WAT"])
SelectorMol( size=630
0: Molecule( WAT:3 num_atoms=3 num_residues=1 )
1: Molecule( WAT:4 num_atoms=3 num_residues=1 )
2: Molecule( WAT:5 num_atoms=3 num_residues=1 )
3: Molecule( WAT:6 num_atoms=3 num_residues=1 )
4: Molecule( WAT:7 num_atoms=3 num_residues=1 )
...
625: Molecule( WAT:628 num_atoms=3 num_residues=1 )
626: Molecule( WAT:629 num_atoms=3 num_residues=1 )
627: Molecule( WAT:630 num_atoms=3 num_residues=1 )
628: Molecule( WAT:631 num_atoms=3 num_residues=1 )
629: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
Note
Note how the index operator will convert the result of the search to
the largest possible container (in this case, Molecule), while
the .residues function will always return the result as a
set of Residue objects.
Another route is to use the atoms in, residues in, chains in, or
segments in phrases in the search, e.g.
>>> print(mols["atoms in molname ACE"])
Selector<SireMol::Atom>( size=22
0: Atom( HH31:1 [ 18.45, 3.49, 12.44] )
1: Atom( CH3:2 [ 18.98, 3.45, 13.39] )
2: Atom( HH32:3 [ 20.05, 3.63, 13.29] )
3: Atom( HH33:4 [ 18.80, 2.43, 13.73] )
4: Atom( C:5 [ 18.48, 4.55, 14.35] )
...
17: Atom( H:18 [ 15.34, 5.45, 17.96] )
18: Atom( CH3:19 [ 13.83, 3.94, 18.35] )
19: Atom( HH31:20 [ 14.35, 3.41, 19.15] )
20: Atom( HH32:21 [ 13.19, 4.59, 18.94] )
21: Atom( HH33:22 [ 13.21, 3.33, 17.69] )
)
Note
Note how the index operator will return a single molecule Selector_Atom_
when only a single molecule matches the search.
>>> print(mols["residues in molnum 10:21"])
SelectorMol( size=11
0: Molecule( WAT:10 num_atoms=3 num_residues=1 )
1: Molecule( WAT:11 num_atoms=3 num_residues=1 )
2: Molecule( WAT:12 num_atoms=3 num_residues=1 )
3: Molecule( WAT:13 num_atoms=3 num_residues=1 )
4: Molecule( WAT:14 num_atoms=3 num_residues=1 )
...
6: Molecule( WAT:16 num_atoms=3 num_residues=1 )
7: Molecule( WAT:17 num_atoms=3 num_residues=1 )
8: Molecule( WAT:18 num_atoms=3 num_residues=1 )
9: Molecule( WAT:19 num_atoms=3 num_residues=1 )
10: Molecule( WAT:20 num_atoms=3 num_residues=1 )
)
You can also go back to the containing molecule using molecules with,
>>> print(mols["molecules with element C"])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
>>> print(mols["molecules with resname /wat/i"])
SelectorMol( size=630
0: Molecule( WAT:3 num_atoms=3 num_residues=1 )
1: Molecule( WAT:4 num_atoms=3 num_residues=1 )
2: Molecule( WAT:5 num_atoms=3 num_residues=1 )
3: Molecule( WAT:6 num_atoms=3 num_residues=1 )
4: Molecule( WAT:7 num_atoms=3 num_residues=1 )
...
625: Molecule( WAT:628 num_atoms=3 num_residues=1 )
626: Molecule( WAT:629 num_atoms=3 num_residues=1 )
627: Molecule( WAT:630 num_atoms=3 num_residues=1 )
628: Molecule( WAT:631 num_atoms=3 num_residues=1 )
629: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
This last search is particularly useful if you are looking for a protein
or for a ligand molecule in a system. In these cases, you could look for
molecules with resname /ala/i or molecules with resname /lig/i.
Uniquely identifying a molecule¶
Molecules are uniquely identified by their molecule number. This is a number
that sire assigns to molecules when they are loaded. sire ensures
that each molecule that it loads will have its own unique number. You have
no control over the number, and should not assume that the number will be
the same every time you run your script. The number is the count of the
molecules that have been loaded during a sire session, e.g. molecule
number 2 refers to the second molecule that has been loaded. You can assume
that the first molecule loaded from a file will have the smallest molecule
number.
You can get all of the numbers of the molecules in a container using
>>> print(mols.numbers())
[MolNum(2), MolNum(3), MolNum(4), MolNum(5)... MolNum(632)]
Note
You can use mols.names() to get all of the molecule names.
sire uses the molecule number as it is the only identifier that can
be guaranteed to be unique in the program. Molecules can be moved and copied
into multiple containers, and its index will depend on its container.
Molecule names are assigned from the file, and multiple molecules can
have the same name.
Molecule identifying types¶
Another way to index molecules is to use the molecule identifying types, i.e.
MolName, MolNum and
MolIdx. The easiest way to create these is via the
function sire.molid().
>>> print(mols[sr.molid("ACE")])
Molecule( ACE:2 num_atoms=22 num_residues=3 )
This returns the molecule called “ACE”, as sr.molid("ACE") has created
an MolName object.
>>> print(sr.molid("ACE"))
MolName('ACE')
This function will create an MolNum if it is passed
an integer, e.g.
>>> print(sr.molid(5))
MolNum(5)
>>> print(mols[sr.molid(5)])
Molecule( WAT:5 num_atoms=3 num_residues=1 )
You can set both a name and a number by passing in two arguments, e.g.
>>> print(mols[sr.molid("WAT", 5)])
Molecule( WAT:5 num_atoms=3 num_residues=1 )
Iterating over molecules¶
The SelectorMol and System classes
are iterable, meaning that they can be used in loops.
>>> for mol in mols["molname WAT and molnum < 10"]:
... print(mol)
Molecule( WAT:3 num_atoms=3 num_residues=1 )
Molecule( WAT:4 num_atoms=3 num_residues=1 )
Molecule( WAT:5 num_atoms=3 num_residues=1 )
Molecule( WAT:6 num_atoms=3 num_residues=1 )
Molecule( WAT:7 num_atoms=3 num_residues=1 )
Molecule( WAT:8 num_atoms=3 num_residues=1 )
Molecule( WAT:9 num_atoms=3 num_residues=1 )
This is particulary useful when combined with looping over the atoms or residues in the molecules.
>>> for mol in mols["molidx < 3"]:
... for atom in mol["element O"]:
... print(mol, atom)
Molecule( ACE:2 num_atoms=22 num_residues=3 ) Atom( O:6 [ 19.19, 5.44, 14.76] )
Molecule( ACE:2 num_atoms=22 num_residues=3 ) Atom( O:16 [ 14.94, 3.17, 15.88] )
Molecule( WAT:3 num_atoms=3 num_residues=1 ) Atom( O:23 [ 25.64, 8.50, 22.42] )
Molecule( WAT:4 num_atoms=3 num_residues=1 ) Atom( O:26 [ 22.83, 8.93, 4.14] )
Smart search terms¶
There are a few smart search terms that can help find molecules.
The most useful is perhaps water. This searches for molecules that
contain one oxygen atom, two hydrogen atoms and any number of null
element (dummy) atoms.
>>> print(mols["water"])
SelectorMol( size=630
0: Molecule( WAT:3 num_atoms=3 num_residues=1 )
1: Molecule( WAT:4 num_atoms=3 num_residues=1 )
2: Molecule( WAT:5 num_atoms=3 num_residues=1 )
3: Molecule( WAT:6 num_atoms=3 num_residues=1 )
4: Molecule( WAT:7 num_atoms=3 num_residues=1 )
...
625: Molecule( WAT:628 num_atoms=3 num_residues=1 )
626: Molecule( WAT:629 num_atoms=3 num_residues=1 )
627: Molecule( WAT:630 num_atoms=3 num_residues=1 )
628: Molecule( WAT:631 num_atoms=3 num_residues=1 )
629: Molecule( WAT:632 num_atoms=3 num_residues=1 )
)
You can also search for protein molecules using protein, e.g.
>>> mols = sr.load("7SA1")
>>> print(mols["protein"])
Molecule( 7SA1:633 num_atoms=11728 num_residues=1518 )
This searches for all molecules that contain at least five residues that are named one of the standard amino acid names.
You can see which names are matched using;
>>> print(sr.search.get_protein_residue_names())
['cyx', 'gly', 'val', 'thr', 'hip', 'his', 'tyr', 'ile', 'trp',
'ala', 'pro', 'glh', 'ash', 'lys', 'ser', 'gln', 'arg', 'asn',
'asp', 'cys', 'met', 'phe', 'leu', 'glu', 'hid', 'hie']
and you can see the minimum number of residues that should be matched using
>>> print(sr.search.get_min_protein_residues())
5
You can set these parameters using the set_protein_residue_names()
and set_min_protein_residues() functions.
Note
Note that the residue names are case insensitive, i.e. ala will match ala, Ala, ALA etc.