0
$\begingroup$

The task I'm thinking of is to detect entire molecules from structures, such as in a crystallographic information file (CIF) or from the XYZ file generated in a molecular dynamics simulation. Crystallography itself doesn't help identify the molecular identity (such as represented in IUPAC name or as a SMILES string), but requires human interpretation or computational modelling.

Three example use cases:

  1. Identify the spatial orientation of all the moieties of molecule X in a molecular crystal of X.
  2. Identify locations of solvent molecules trapped in crystals
  3. Identify parts of a molecule that has multiple copies in the crystal.

Does anyone know of such an algorithm/software?

$\endgroup$
9
  • 2
    $\begingroup$ If you can get the structure why would you not know the name beforehand? $\endgroup$
    – Mithoron
    Jun 16, 2022 at 15:52
  • $\begingroup$ I'm looking for part of a structure that is still molecular. In other scenarios, if you have multiple copies of the same molecule, their spatial orientation (such as packing) would be interesting to understand. Neither of these would be possible without identifying the molecular entity in 3D (e.g. as a grouping) $\endgroup$
    – QGent
    Jun 16, 2022 at 16:06
  • 1
    $\begingroup$ Its not clear what you want as the cif/sdf/pdb etc. data (obtained usually from x-ray crystallography) contains all the information, i.e. xyz coords and connection table to work out where the bonds are, so you can draw the molecule at a pinch by hand from the connection table. If you don't have a connection table you can use known bond lengths and atom positions to work one out by writing some simple code. $\endgroup$
    – porphyrin
    Jun 17, 2022 at 11:04
  • 1
    $\begingroup$ Sorry, I don’t think that helps much. Perhaps if you have an example CIF in which you think this is necessary? The space group, for example will give you all the symmetric copies and orientations of a molecule in the unit cell. Finding a solvent isn’t very hard either, because you can look through the bonded components to find solvents or co-crystals. $\endgroup$ Jun 18, 2022 at 20:03
  • 1
    $\begingroup$ Load the file into a 3D browser like JSmol. It will show bonds where they make sense, allowing you to distinguish between distinct molecules in the structure. Some cases might be more complex and you would need to know what is considered a molecule in that context. $\endgroup$
    – Karsten
    Jun 18, 2022 at 21:44

1 Answer 1

2
$\begingroup$

Your question has a few parts, but multiple solutions. Let's start with "XYZ from molecular dynamics"

3
XYZ file generated by Avogadro.
O     -2.86063    1.04260   -0.00000
H     -2.52287    1.38117    0.82085
H     -3.53528    0.39764    0.17713

This is trivial. Open Babel will read this, perceive the bonds through inter-atomic distances, and generate water. Most programs do this. (RDKit should use the xyz2mol package.)

Problems come when you want to perceive bond orders .. because it's a non-trivial task. Both Open Babel and xyz2mol have algorithms..

You ask about solvent molecules. Again, both Open Babel and RDKit have methods to find and/or remove fragments, e.g. "An open source chemical structure curation pipeline using RDKit"

Then you ask about multiple copies (e.g., in a CIF)

Here's part of an aspirin entry in the COD:

O1 O 0.62354(11) 0.1419(2) 0.61266(11) 0.0203(3) Uani 1 1 d .
H1 H 0.577(2) 0.025(5) 0.602(2) 0.060(8) Uiso 1 1 d .
O2 O 0.50970(11) 0.1879(2) 0.40933(11) 0.0212(3) Uani 1 1 d .
O3 O 0.78878(10) 0.41876(19) 0.72945(10) 0.0171(3) Uani 1 1 d .
O4 O 0.90669(11) 0.2189(2) 0.66325(11) 0.0215(3) Uani 1 1 d .
C1 C 0.65418(15) 0.4419(3) 0.50750(16) 0.0154(4) Uani 1 1 d .
C2 C 0.74806(15) 0.5203(3) 0.61181(16) 0.0162(4) Uani 1 1 d .
C3 C 0.80069(16) 0.7061(3) 0.60519(17) 0.0196(4) Uani 1 1 d .
H3 H 0.8632 0.7580 0.6776 0.024 Uiso 1 1 calc R
C4 C 0.76267(16) 0.8174(3) 0.49303(18) 0.0215(4) Uani 1 1 d .
H4 H 0.7997 0.9449 0.4882 0.026 Uiso 1 1 calc R
C5 C 0.67053(16) 0.7425(3) 0.38781(17) 0.0195(4) Uani 1 1 d .
H5 H 0.6444 0.8182 0.3106 0.023 Uiso 1 1 calc R
C6 C 0.61699(15) 0.5578(3) 0.39559(16) 0.0173(4) Uani 1 1 d .
H6 H 0.5534 0.5081 0.3233 0.021 Uiso 1 1 calc R
C7 C 0.59016(15) 0.2458(3) 0.50601(15) 0.0151(4) Uani 1 1 d .
C8 C 0.86823(15) 0.2623(3) 0.74362(16) 0.0166(4) Uani 1 1 d .
C9 C 0.89840(17) 0.1611(3) 0.86922(17) 0.0249(5) Uani 1 1 d .
H7 H 0.9545 0.0480 0.8770 0.037 Uiso 1 1 calc R
H8 H 0.8259 0.1059 0.8769 0.037 Uiso 1 1 calc R
H9 H 0.9347 0.2625 0.9366 0.037 Uiso 1 1 calc R

That's it. That's the connected component for an aspirin molecule. If there were solvent, it would be a disconnected fragment if you, e.g. convert to SDF or SMILES.

How do you know the symmetric copies? It's in the CIF as defined by the space group.

_symmetry_equiv_pos_as_xyz
'x, y, z'
'-x, y+1/2, -z+1/2'
'-x, -y, -z'
'x, -y-1/2, z-1/2'

This gives you the transforms (in fractional coordinates) for all the symmetric copies. Some may be duplicates. But if you have a point $(x,y,z)$ there's a symmetric copy at $(-x, y + 0.5, -z + 0.5)$, etc. In other words this space group could have 4 copies of aspirin per unit cell.

filled unit cell of aspirin

You mention automation, which sounds like you want to write a script. Maybe that's in Open Babel. Maybe you use CCDC Python or ASE or pymatgen. There are likely other packages for handling CIF files. I don't know what precise features you want. But they should all find connected components. Once you've built the script to process one, adding an outer loop in Python to handle multiple files is fairly easy, e.g.:

for filename in glob.iglob("*.cif"):
    # do something
$\endgroup$
1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.