# How to sort a supercell for a molecular crystal?

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_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_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_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

• what does sort by molecule mean? Mar 30, 2021 at 1:58
• Sadly, if you don't answer Cody's question, your \$100 bounty may go wasted. Why not do the most you can to try to get the best value out of your \$100 bounty? Apr 1, 2021 at 20:54
• @CodyAldaz By "sort by molecule," I mean all of the atoms for molecule 1, in a specified order, followed by all of the atoms in the same specified order, for molecule 2, and so on. Let me know if that's not clear. Apr 1, 2021 at 22:57

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

The following code does this:

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

# 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

• +1. nice! I hope this solves the OPs problem since they spend a \$100 bounty on it! By the way, did you know we have an ASE chat room? It would be nice to have more people participating there too! Apr 2, 2021 at 19:45
• Wonderful! Thank you so much! I'm surprised that it isn't easier in all of these software packages to do this. Apr 2, 2021 at 22:21

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!

• +1 Thanks for taking care of this one Matt! By the way, one of are longest standing unanswered questions is about pymatgen, do you have any idea how to answer it? Mar 12, 2021 at 20:22
• @MattHorton This looks good, but how do you sort a supercell that is generated from the unit cell? And how would you reverse the order given by the last line? From what I can tell, you only sorted one molecule. Mar 13, 2021 at 15:05
• How do I do this for a supercell? I tried including acetonitrile = acetonitrile * (5,5,5) but it prints only one molecule when I type print(molecules), and I also don't know how to sort all of the molecules in the supercell in the reverse order that you did for a single molecule. Mar 14, 2021 at 23:28
• @user1, in pymatgen there is a Structure object (for periodic crystals) and a Molecule object (for isolated molecules), the Structure can be sorted in a similar way with its own from_sites method. The most difficult thing in your case is figuring out what logic you want to use for the sorting. Mar 15, 2021 at 19:44
• @NikeDattani, do you have a link? Happy to answer any pymatgen questions, they're more likely to get a rapid response via our own support forum however. Mar 15, 2021 at 19:45