How to sort a supercell for a molecular crystal?

Question

I am trying to build a supercell of the following (below) molecular crystal from its .cif file. How do I sort a supercell here according to each molecule in the supercell? I'm not tethered to a specific software to do this, just something that works as I want it to.

I have tried using ASE, Vesta, and Avogadro to do this, without success, along with the pymatgen suggestion below. It seems that these software are better set up to do this for an ionic crystal rather than a supercell of a molecule. For example, in the Atomic Simulation Environment (ASE), when I use get_chemical_symbols(), I see more atoms than there are in a single molecule. Furthermore, I do not know how to distinguish between the carbon bonded to the hydrogens and the central carbon (molecule is CH3CN). What I want is an output file (say, .xyz) with all of the atoms, in a specified order (say, H, H, H, methyl carbon, nitrile carbon, N), for molecule 1 in the supercell, followed by all of the atoms in the same specified order, for molecule 2, and so on.

data_I
_audit_creation_method     SHELXL
_journal_date_recd_electronic     2002-06-07
_journal_date_accepted     2002-09-25
_journal_name_full     'Acta Crystallographica, Section B'
_journal_year     2002
_journal_volume     58
_journal_issue      6
_journal_page_first     1005
_journal_page_last     1010
_journal_paper_category     FA
_chemical_name_systematic     'acetonitrile'
_chemical_name_common     ?
_chemical_formula_moiety     'C2 H3 N'
_chemical_formula_sum     'C2 H3 N'
_chemical_formula_structural     'C H3 C N'
_chemical_formula_analytical     ?
_chemical_formula_weight     41.05
_chemical_melting_point     ?
_symmetry_cell_setting     'monoclinic'
_symmetry_space_group_name_H-M     'P 21/c'
_symmetry_space_group_name_Hall     '-P 2ybc'
loop_
    _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'
_cell_length_a     4.102(3)
_cell_length_b     8.244(7)
_cell_length_c     7.970(7)
_cell_angle_alpha     90.00
_cell_angle_beta     100.10(10)
_cell_angle_gamma     90.00
_cell_volume     265.3(4)
_cell_formula_units_Z     4
_cell_measurement_reflns_used     25
_cell_measurement_theta_min     5.20
_cell_measurement_theta_max     21.48
_cell_measurement_temperature     201(2)
_exptl_crystal_description     'cylinder'
_exptl_crystal_colour     'colourless'
_exptl_crystal_size_max     1.2
_exptl_crystal_size_mid     0.5
_exptl_crystal_size_min     0.3
_exptl_crystal_size_rad     0.15
_exptl_crystal_density_diffrn     1.028
_exptl_crystal_density_meas     ?
_exptl_crystal_density_method     ?
_exptl_crystal_F_000     88
_exptl_absorpt_coefficient_mu     0.067
_exptl_absorpt_correction_type     'none'
_exptl_absorpt_correction_T_min     ?
_exptl_absorpt_correction_T_max     ?
_exptl_special_details
;
 ?
;
_diffrn_ambient_temperature     201(2)
_diffrn_radiation_type     MoK\a
_diffrn_radiation_wavelength     0.71073
_diffrn_radiation_source     'fine-focus sealed tube'
_diffrn_radiation_monochromator     'graphite'
_diffrn_measurement_device     'Nonius CAD4 diffractometer'
_diffrn_measurement_method     '\w--2\q'
_diffrn_reflns_number     376
_diffrn_reflns_av_R_equivalents     0.0571
_diffrn_reflns_av_sigmaI/netI     0.0591
_diffrn_reflns_theta_min     3.58
_diffrn_reflns_theta_max     21.89
_diffrn_reflns_theta_full     21.89
_diffrn_measured_fraction_theta_max     0.99
_diffrn_measured_fraction_theta_full     0.99
_diffrn_reflns_limit_h_min     0
_diffrn_reflns_limit_h_max     4
_diffrn_reflns_limit_k_min     0
_diffrn_reflns_limit_k_max     8
_diffrn_reflns_limit_l_min     -8
_diffrn_reflns_limit_l_max     8
_diffrn_standards_number     'none'
_diffrn_standards_interval_count     'none'
_diffrn_standards_interval_time     'none'
_diffrn_standards_decay_%     'none'
_refine_special_details
;
 Refinement of F^2^ against ALL reflections.  Weighted R-factors wR and
 goodnesses of fit S are based on F^2^, conventional R-factors R are based
 on F, with F set to zero for negative F^2^. The threshold_expression of
 F^2^ > 2sigma(F^2^) is used only for calculating R_factors(gt) etc. and is
 not relevant to the choice of reflections for refinement.  R-factors based
 on F^2^ are statistically about twice as large as those based on F, and R-
 factors based on ALL data will be even larger.
;
_reflns_number_total     324
_reflns_number_gt     202
_reflns_threshold_expression     'I>2\s(I)'
_refine_ls_structure_factor_coef     Fsqd
_refine_ls_matrix_type     full
_refine_ls_R_factor_all     0.1051
_refine_ls_R_factor_gt     0.0472
_refine_ls_wR_factor_all     0.1504
_refine_ls_wR_factor_ref     0.1106
_refine_ls_goodness_of_fit_all     1.137
_refine_ls_goodness_of_fit_ref     1.107
_refine_ls_restrained_S_all     1.137
_refine_ls_restrained_S_obs     1.107
_refine_ls_number_reflns     324
_refine_ls_number_parameters     41
_refine_ls_number_restraints     0
_refine_ls_hydrogen_treatment     'refall'
_refine_ls_weighting_scheme     calc
_refine_ls_weighting_details
              'w=1/[\s^2^(Fo^2^)+(0.0360P)^2^+0.1879P] where P=(Fo^2^+2Fc^2^)/3'
_atom_sites_solution_hydrogens     difmap
_atom_sites_solution_primary     direct
_atom_sites_solution_secondary     difmap
_refine_ls_shift/su_max     0.000
_refine_ls_shift/su_mean     0.000
_refine_diff_density_max     0.169
_refine_diff_density_min     -0.168
_refine_ls_extinction_method     SHELXL
_refine_ls_extinction_coef     0.07(4)
_refine_ls_extinction_expression
                                  'Fc^*^=kFc[1+0.001xFc^2^\l^3^/sin(2\q)]^-1/4^'
loop_
    _atom_type_symbol
    _atom_type_description
    _atom_type_scat_dispersion_real
    _atom_type_scat_dispersion_imag
    _atom_type_scat_source
    'C' 'C' 0.0033 0.0016
                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    'H' 'H' 0.0000 0.0000
                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
    'N' 'N' 0.0061 0.0033
                         'International Tables Vol C Tables 4.2.6.8 and 6.1.1.4'
_computing_data_collection     'CAD4-EXPRESS (Enraf-Nonius, 1993)'
_computing_cell_refinement     'CAD4-EXPRESS (Enraf-Nonius, 1993)'
_computing_data_reduction     'CADAK (Savariault,1991)'
_computing_structure_solution     'SHELXS-96 (Sheldrick, 1990)'
_computing_structure_refinement     'SHELXL-96 (Sheldrick, 1996)'
_computing_molecular_graphics     'ORTEP III (Burnett & Johnson, 1996)'
_computing_publication_material     'SHELXL-96 (Sheldrick, 1996)'
loop_
    _atom_site_label
    _atom_site_fract_x
    _atom_site_fract_y
    _atom_site_fract_z
    _atom_site_U_iso_or_equiv
    _atom_site_thermal_displace_type
    _atom_site_calc_flag
    _atom_site_refinement_flags
    _atom_site_occupancy
    _atom_site_disorder_group
    _atom_site_type_symbol
    N 0.4547(9) 0.2657(5) 0.4613(4) 0.0710(15) Uani d . 1 . N
    C1 0.0949(12) 0.4579(6) 0.2478(6) 0.0586(14) Uani d . 1 . C
    C2 0.2946(9) 0.3498(5) 0.3672(5) 0.0507(13) Uani d . 1 . C
    H1 -0.108(11) 0.402(5) 0.166(5) 0.089(14) Uiso d . 1 . H
    H2 0.233(11) 0.518(6) 0.186(6) 0.113(18) Uiso d . 1 . H
    H3 -0.050(11) 0.538(6) 0.301(5) 0.103(16) Uiso d . 1 . H
loop_
    _atom_site_aniso_label
    _atom_site_aniso_U_11
    _atom_site_aniso_U_22
    _atom_site_aniso_U_33
    _atom_site_aniso_U_12
    _atom_site_aniso_U_13
    _atom_site_aniso_U_23
    N 0.076(3) 0.070(3) 0.064(2) 0.005(2) 0.0028(18) 0.008(2)
    C1 0.058(3) 0.060(3) 0.055(3) 0.005(2) 0.001(2) 0.008(2)
    C2 0.056(2) 0.050(3) 0.047(2) -0.007(2) 0.0101(19) -0.009(2)
_geom_special_details
;
 All esds (except the esd in the dihedral angle between two l.s. planes)
 are estimated using the full covariance matrix.  The cell esds are taken
 into account individually in the estimation of esds in distances, angles
 and torsion angles; correlations between esds in cell parameters are only
 used when they are defined by crystal symmetry.  An approximate (isotropic)
 treatment of cell esds is used for estimating esds involving l.s. planes.
;
loop_
    _geom_bond_atom_site_label_1
    _geom_bond_atom_site_label_2
    _geom_bond_site_symmetry_2
    _geom_bond_distance
    _geom_bond_publ_flag
    N C2 . 1.141(5) yes
    C1 C2 . 1.448(6) yes
    C1 H1 . 1.07(5) yes
    C1 H2 . 0.96(5) yes
    C1 H3 . 1.03(5) yes
loop_
    _geom_angle_atom_site_label_1
    _geom_angle_atom_site_label_2
    _geom_angle_atom_site_label_3
    _geom_angle_site_symmetry_1
    _geom_angle_site_symmetry_3
    _geom_angle
    _geom_angle_publ_flag
    C2 C1 H1 . . 115(2) yes
    C2 C1 H2 . . 110(3) yes
    H1 C1 H2 . . 112(3) yes
    C2 C1 H3 . . 115(2) yes
    H1 C1 H3 . . 95(3) yes
    H2 C1 H3 . . 108(4) yes
    N C2 C1 . . 179.3(4) yes

Answer 1

ASE can be used to achieve what you require. Basically we need to:

• Read the structure.
• Determine the molecules using neighborlist class in ASE.
• Sort each of the molecules in the required manner using ase.build.sort.

The following code does this:

import numpy as np
from ase.io import read, write
from ase import Atoms, neighborlist
from ase.build import sort
from scipy import sparse                                                   

acetonitrile = read("acetonitrile.cif") # Read the structure

# Determine the indices of each molecule using neighborlists
cutoff = neighborlist.natural_cutoffs(acetonitrile)
nl = neighborlist.build_neighbor_list(acetonitrile, cutoffs=cutoff)
connmat = nl.get_connectivity_matrix(False) # Connectivity matrix
# n_components contains number of molecules
# component_list contains molecule number of each atom in the system
n_components, component_list = sparse.csgraph.connected_components(connmat)

# Sorting based on the molecule index first (`component_list`)
# and then by chemical symbol of the atom
sortedacetonitrile = sort(
        acetonitrile,
        tags=[(component_list[i]*10) + ['H','C','N'].index(s)
            for i,s in enumerate(acetonitrile.get_chemical_symbols())]
        )

# Write to file
write("sortedacetonitrile.xyz", sortedacetonitrile)

Code to determine molecules was taken from an example in the ASE website.

The resulting XYZ file:

24
Lattice="4.102 0.0 0.0 0.0 8.244 0.0 -1.3976728069531361 0.0 7.84649034439626" Properties=species:S:1:pos:R:3:spacegroup_kinds:I:1 spacegroup="P 21/c" unit_cell=conventional occupancy="_JSON {\"0\": {\"N\": 1}, \"1\": {\"C\": 1}, \"2\": {\"C\": 1}, \"3\": {\"H\": 1}, \"4\": {\"H\": 1}, \"5\": {\"H\": 1}}" pbc="T T T"
H        3.42697031       3.31408800       1.30251740        3
H        0.69579886       4.27039200       1.45944720        4
H        3.47620049       4.43527200       2.36179359        5
C        0.04293648       3.77492760       1.94436031        1
C        0.69522375       2.88375120       2.88123125        2
N        1.22043293       2.19043080       3.61958600        0
H       -0.02380672       7.43608800       2.62072778        3
H        2.70736474       0.14839200       2.46379797        4
H       -0.07303689       0.31327200       1.56145158        5
C        3.36022712       7.89692760       1.97888486        1
C        2.70793985       7.00575120       1.04201392        2
N        2.18273066       6.31243080       0.30365918        0
H       -0.72264312       4.92991200       6.54397295        3
H        2.00852834       3.97360800       6.38704314        4
H       -0.77187329       3.80872800       5.48469675        5
C        2.66139071       4.46907240       5.90213004        1
C        2.00910345       5.36024880       4.96525909        2
N        1.48389426       6.05356920       4.22690435        0
H        2.72813391       0.80791200       5.22576257        3
H       -0.00303755       8.09560800       5.38269238        4
H        2.77736408       7.93072800       6.28503877        5
C       -0.65589993       0.34707240       5.86760548        1
C       -0.00361266       1.23824880       6.80447643        2
N        0.52159653       1.93156920       7.54283117        0

Answer 2

I don't know if this is exactly what you want to do, but pymatgen may be able to help (disclosure, I am a developer of this package).

See this example code below for your example input file:

from pymatgen.core.structure import Structure, Molecule
from pymatgen.analysis.graphs import StructureGraph
from pymatgen.analysis.local_env import JmolNN

# load the molecular crystal from a CIF file
# if a single CIF file contains multiple structures, 
# you will need to use pymatgen.io.cif.CifParser(your_file).get_structures() instead
acetonitrile = Structure.from_file("acetonitrile.cif")

# attempt to guess the bonds present using a bonding strategy
# (here JmolNN but others are available)
# this creates a StructureGraph
sg = StructureGraph.with_local_env_strategy(acetonitrile, JmolNN())

# now extract individual molecules, which are isolated subgraphs
# (i.e. groups of atoms not connected to each other)
molecules = sg.get_subgraphs_as_molecules()

# this gives a list of molecules
# here, there is only a single unique subgraph so only one is returned

print(molecules[0])
# Molecule Summary
# Site: H (-1.5303, 0.3235, -1.2207)
# Site: C (-0.7109, 0.7843, -0.7147)
# Site: H (-1.2961, 1.4446, -0.1864)
# Site: C (0.0955, -0.1069, 0.0933)
# Site: N (0.7421, -0.8002, 0.7281)
# Site: H (-0.1532, 1.2798, -1.3066)

# you can then sort the sites according to whatever logic you want,
# and create a new Molecule object from the sorted sites
new_molecule = Molecule.from_sites(sorted(molecules[0], key=lambda site: site.species))
print(new_molecule)

# Molecule Summary
# Site: N (0.7421, -0.8002, 0.7281)
# Site: C (-0.7109, 0.7843, -0.7147)
# Site: C (0.0955, -0.1069, 0.0933)
# Site: H (-1.5303, 0.3235, -1.2207)
# Site: H (-1.2961, 1.4446, -0.1864)
# Site: H (-0.1532, 1.2798, -1.3066)

Hope this helps!