8
$\begingroup$

If I have an XYZ file with a cluster of a dozen of molecules: water, $\ce{O2}$, $\ce{H2}$, $\ce{H2O2}$. A human may easily identify these molecules. But how to do it automatically? I need to identify let's say all water and $\ce{H2}$ molecules, and I know nothing about other molecules ($\ce{O2}$, $\ce{H2O2}$). Generally speaking, I need to recognize a few known small molecules in arbitrary systems of medium size (hundreds of atoms). These systems may contain unknown molecules, or even bulk materials. Please recommend me some solution for that task. I would be happy to find either some existing software (like Python package), or some not-too-complex algorithm that I could implement myself.

Here is an exemplar structure in XYZ format:

20

O     -1.47655    1.21497   -0.20250
H     -0.53716    1.39872   -0.20387
H     -1.84821    1.86921    0.38917
O     -1.20797   -1.49263   -0.10100
H     -1.40154   -0.55535   -0.11741
H     -1.71280   -1.85472   -0.82920
O      1.45295   -1.19622    0.39415
H      0.53738   -1.39890    0.20209
H      1.65598   -1.71658    1.17148
O      1.23158    1.47389   -0.09064
H      1.40132    0.55553    0.11919
H      1.90503    1.70209   -0.73144
O     -0.76633    3.57210   -0.37181
O      0.15413    3.64478   -0.39400
O      2.01459    3.79209   -0.43890
O      3.44119    3.98651   -0.48207
H      3.86223    3.37474    0.11035
H      1.59355    4.40386   -1.03132
H     -4.39535    2.45856   -0.19546
H     -4.82240    2.33535   -0.65595
$\endgroup$

5 Answers 5

9
$\begingroup$

I used OpenBabel with the --separate option. From the documentation:

--separate
Separate disconnected fragments into individual molecular records

I created a file called water_cluster.xyz:

20

O     -1.47655    1.21497   -0.20250
H     -0.53716    1.39872   -0.20387
H     -1.84821    1.86921    0.38917
O     -1.20797   -1.49263   -0.10100
H     -1.40154   -0.55535   -0.11741
H     -1.71280   -1.85472   -0.82920
O      1.45295   -1.19622    0.39415
H      0.53738   -1.39890    0.20209
H      1.65598   -1.71658    1.17148
O      1.23158    1.47389   -0.09064
H      1.40132    0.55553    0.11919
H      1.90503    1.70209   -0.73144
O     -0.76633    3.57210   -0.37181
O      0.15413    3.64478   -0.39400
O      2.01459    3.79209   -0.43890
O      3.44119    3.98651   -0.48207
H      3.86223    3.37474    0.11035
H      1.59355    4.40386   -1.03132
H     -4.39535    2.45856   -0.19546
H     -4.82240    2.33535   -0.65595

and run the command:

obabel -iXYZ water_cluster.xyz --separate -oXYZ

The output was:

3
water_cluster.xyz#1
O         -1.47655        1.21497       -0.20250
H         -0.53716        1.39872       -0.20387
H         -1.84821        1.86921        0.38917
3
water_cluster.xyz#2
O         -1.20797       -1.49263       -0.10100
H         -1.40154       -0.55535       -0.11741
H         -1.71280       -1.85472       -0.82920
3
water_cluster.xyz#3
O          1.45295       -1.19622        0.39415
H          0.53738       -1.39890        0.20209
H          1.65598       -1.71658        1.17148
3
water_cluster.xyz#4
O          1.23158        1.47389       -0.09064
H          1.40132        0.55553        0.11919
H          1.90503        1.70209       -0.73144
2
water_cluster.xyz#5
O         -0.76633        3.57210       -0.37181
O          0.15413        3.64478       -0.39400
4
water_cluster.xyz#6
O          2.01459        3.79209       -0.43890
O          3.44119        3.98651       -0.48207
H          3.86223        3.37474        0.11035
H          1.59355        4.40386       -1.03132
2
water_cluster.xyz#7
H         -4.39535        2.45856       -0.19546
H         -4.82240        2.33535       -0.65595

Adding -m to the command above, and an output file name:

obabel -iXYZ water_cluster.xyz --separate -m -oXYZ -O fragment.xyz

will produce separated files, one for each fragment:

-a----        2024-04-04     16:15            169 fragment1.xyz
-a----        2024-04-04     16:15            169 fragment2.xyz
-a----        2024-04-04     16:15            169 fragment3.xyz
-a----        2024-04-04     16:15            169 fragment4.xyz
-a----        2024-04-04     16:15            120 fragment5.xyz
-a----        2024-04-04     16:15            218 fragment6.xyz
-a----        2024-04-04     16:15            120 fragment7.xyz
$\endgroup$
1
  • $\begingroup$ How does it work? Does it identify the molecules? I need to couple these fragments with known molecules. In this case it can be done via stoichiometry, but how about more complex cases? Can it be used with ASE.Atoms ? $\endgroup$
    – user36313
    Apr 5 at 10:12
7
$\begingroup$

You ask for either software or algorithm. Other answers speak to existing software. The underlying algorithms generally revolve around interatomic distances.

For instance,

  • one chooses a suitable set of atomic radii for the elements involved, and

  • defines a bond to exist between each pair of atoms that are closer than the sum of their radii, plus a smallish fudge factor to account for stretched bonds and imprecise positions.

  • possibly one then prunes the connectivity list in any of a variety of ways, such as by removing all but the shortest of each hydrogen atom's bonds (to split hydrogen-bonded pairs of molecules) or similar.

  • Having thereby identified the molecular connectivity, one can fairly easily split the connectivity graph into its maximal connected subgraphs. Each of these corresponds to a molecule, an ion, or perhaps to a polymer or a metalic or ionic cluster.

One can even estimate bond order from the results, as this is pretty strongly correlated to bond length, at least for first-row elements.

$\endgroup$
3
$\begingroup$

You can use some clustering algorithm, or a mixture model.

E.g., in Julia:

using Clustering

This will import the Clustering.jl library.

Then, a function to read your data:

function read_xyz(input_xyz)
    lines = open(input_xyz) do file; readlines(file); end

    # Number of atoms
    n_atoms = parse(Int, filter(x -> '0' <= x <= '9', lines[1]))

    # Array of atom symbols and coordinates
    atom_symbols = Array{String}(undef, n_atoms)
    coords = Array{Float64}(undef, (3, n_atoms))
    
    # Indices of the first and the last lines with atomic coordinates
    ib = 3; ie = 3 + (n_atoms - 1);
    
    for (i,iline) ∈ enumerate(ib:ie)
        atom_symbols[i] = string(split(lines[iline])[1])
        coords[:,i] .= parse.(Float64, string.(split(lines[iline])[2:4]))
    end
    return atom_symbols, coords
end;

Then, implement a clustering algorithm, e.g., DBSCAN:

let
    input_xyz = "water_cluster.xyz"
    atom_symbol, coords = read_xyz(input_xyz)

    r_db = dbscan(coords, 1.5);

    for i ∈ 1:length(r_db.counts)
        println("Molecule ", i)
        idxs = findall(x -> x == i, r_db.assignments)
        for idx ∈ idxs
            println(atom_symbol[idx], " ", transpose(coords[:,idx]))
        end
        println()
    end 
end

This will print

Molecule 1
O [-1.47655 1.21497 -0.2025]
H [-0.53716 1.39872 -0.20387]
H [-1.84821 1.86921 0.38917]

Molecule 2
O [-1.20797 -1.49263 -0.101]
H [-1.40154 -0.55535 -0.11741]
H [-1.7128 -1.85472 -0.8292]

Molecule 3
O [1.45295 -1.19622 0.39415]
H [0.53738 -1.3989 0.20209]
H [1.65598 -1.71658 1.17148]

Molecule 4
O [1.23158 1.47389 -0.09064]
H [1.40132 0.55553 0.11919]
H [1.90503 1.70209 -0.73144]

Molecule 5
O [-0.76633 3.5721 -0.37181]
O [0.15413 3.64478 -0.394]

Molecule 6
O [2.01459 3.79209 -0.4389]
O [3.44119 3.98651 -0.48207]
H [3.86223 3.37474 0.11035]
H [1.59355 4.40386 -1.03132]

Molecule 7
H [-4.39535 2.45856 -0.19546]
H [-4.8224 2.33535 -0.65595]

See also clustering in Python, R, Java/Scala.

$\endgroup$
2
$\begingroup$

For the task I recommend OpenBabel and a subsequent query on PubChem with PubChemPy. The "dictionary of chemical compounds" used hence is the one by NIH presuming compounds of your cluster are indexed by this database with a matching SMILES string.

In the first step, one wants to identify which atoms are close enough to each other to be considered part of a molecule. Thus, a temporary .sdf file about a "cluster molecule" is split into individual, separate molecules each expressed by its own SMILES string. The following script for the bash shell can be used to process e.g., water_cluster.xyz as input file:

#!/usr/bin/bash

# name: xyz2smi.sh
# use : bash ./xyz2smi.sh input.xyz > output.smi

obabel "$1" -O temp.sdf && \
  obabel temp.sdf -osmi --separate && \
  rm temp.sdf

The intermediate .sdf file considers bond orders different from one

enter image description here

The listing of SMILES in the particular case provided is

O   water_cluster.xyz#1
O   water_cluster.xyz#2
O   water_cluster.xyz#3
O   water_cluster.xyz#4
O=O water_cluster.xyz#5
OO  water_cluster.xyz#6
[H][H]  water_cluster.xyz#7

To identify the chemical names and eventually tally the molecules, PubChem's CACTUS servers are accessed. (Interestingly, PubChem's first chemical name for water is oxidane.)

occurrences     | motif
4       | oxidane
1       | molecular oxygen
1       | hydrogen peroxide
1       | molecular hydrogen
7       | molecules in total

The Python script for this part relies on the non-standard PubChemPy library initiated by Matt Swain equally available via PyPI.

#!/usr/bin/env python3
"""
purpose     : Contact NIH Cactus servers to identify molecules expressed as SMILES string.
inspired by : Oliver Scott, 2020-10-13
source      : https://stackoverflow.com/questions/64329049/converting-smiles-to-chemical-name-or-iupac-name-using-rdkit-or-other-python-mod
"""

import argparse

import pubchempy  # non-standard library initiated by Matt Swain, tested: 1.0.4


def get_args():
    """Get command-line arguments"""

    parser = argparse.ArgumentParser(
        description="""
    Identify and sum up SMILES expressed molecules with NIH Cactus servers""",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )

    parser.add_argument(
        "file",
        help="A list of SMILES strings",
        metavar="FILE",
        type=argparse.FileType("rt", encoding="utf8"),
        default=None,
    )

    return parser.parse_args()


def main():
    """Join the functionalities"""
    args = get_args()
    file_arg = args.file

    compounds_list = []
    for line in file_arg:
        try:
            smiles = line.split()[0]

            compounds = pubchempy.get_compounds(smiles, namespace="smiles")
            match = compounds[0]
            compound = match.iupac_name
            # print(compound)

            compounds_list.append(compound)
        except Exception:
            print(f"failed entry: {smiles}")

    count = {}
    for compound in compounds_list:
        count.setdefault(compound, 0)
        count[compound] = count[compound] + 1

    print("occurrences \t| motif")
    for key, value in count.items():
        print(f"{value} \t\t| {key}")
    print(f"{sum(count.values())} \t\t| molecules in total")


if __name__ == "__main__":
    main()


My first suggestion relied on Jan Jensen's xyz2mol which however for the present cluster data was successful to recover molecular oxygen with a bond order of two, and molecular hydrogen with a bond order of one. Despite that I still think it can be useful for the "recovery" of bonds after a quantum chemical computation (equally see Greg Landrum's corresponding RDKit blog post):

from rdkit import Chem
from rdkit.Chem import rdDetermineBonds

raw_mol = Chem.MolFromXYZFile("water_cluster.xyz")
mol = Chem.Mol(raw_mol)
rdDetermineBonds.DetermineBonds(mol,charge=0)

print(Chem.MolToMolBlock(mol))
$\endgroup$
4
  • 1
    $\begingroup$ Good answer. But i think OP is asking for splitting the file into multiple fragments. $\endgroup$ Apr 5 at 2:10
  • $\begingroup$ And how to identify known molecules here? $\endgroup$
    – user36313
    Apr 5 at 10:13
  • $\begingroup$ @HemanthHaridas The initial answer was edited to account for this point. $\endgroup$
    – robert
    Apr 5 at 22:17
  • $\begingroup$ @user36313 The revised answer now includes a script to seek in NIH Cactus/PubChem servers for chemical names of the molecules identified. Please have a look. $\endgroup$
    – robert
    Apr 5 at 22:18
1
$\begingroup$

You can use the RDKit package to parse your molecule and then iteratively check if a bond is present between each of the atoms. This would allow you to build a connectivity matrix, that you can store as an undirected graph. You can then use the routine strong_components = list(nx.weakly_connected_components(_totalGraph)) to give you all the subgraphs (which in your case will be individual molecules). If you can provide a test xyz file, I can modify the answer to include a python code that can do this.

$\endgroup$
1
  • $\begingroup$ Thanks, I have added the exemplar structure to the OP $\endgroup$
    – user36313
    Apr 4 at 18:45

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .