Are there tools or algorithms to detect organic molecules from 3D crystal structures?

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?




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](https://github.com/jensengroup/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"](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-020-00456-1)


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


Here's part of an [aspirin entry in the COD](http://crystallography.net/cod/1515581.html):



```
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](https://i.stack.imgur.com/3HzaX.png)](https://i.stack.imgur.com/3HzaX.png)


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

```

 2