How to discard molecules at the boundary in the Atomic Simulation Environment (ASE)?

Question

When building a supercell for a molecular crystal with the Atomic Simulation Environment, it seems that the program takes into account periodic boundary conditions when replicating the unit cell. How do I discard all molecules that wrap around the boundaries? (It seems that taking into account the boundary conditions for non-orthorhombic cells is non-trivial, so I'd like to avoid it if possible!)

Please note the question here that this follows up on. The .cif file that I am using is reproduced below.

EDIT: When I try as a test the following:

import ase
from ase.io import read, write
unitcell = read("mycif.cif")
unitcell.set_pbc((False,False,False))
supercell = unitcell*(8,4,4)

I still get unconnected atoms at the boundaries. Any idea about how to resolve this?

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

For an atom in the boundary, the distance with all other atoms of the molecule it is part of will be different based on whether we take into account minimum image convention (MIC) or not. Thus we can,

• first identify the indices of all the molecules in the structure and
• then exclude the molecules for which any of the distance is different if MIC is considered or not.

See the following code for an example of this using ASE. The code to determine indices of molecules is the same as a previous answer.

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

supercell = read("mycif.cif")*(8,4,4) # Read the structure

# Determine the indices of each molecule using neighborlists
cutoff = neighborlist.natural_cutoffs(supercell)
nl = neighborlist.build_neighbor_list(supercell, 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)

# Get symbols array of atoms
origsymbols = np.array(supercell.get_chemical_symbols())

# Excluding atoms in the boundary
ind_selected = [] # List to store selected indices
for i in range(n_components): # For each molecule
    ind = np.where(component_list == i)[0] # Find indices of atoms in molecule
    # Based on MIC, if any atoms are in the boudary exclude the molecule
    if (supercell.get_distances(ind[0],ind) != supercell.get_distances(ind[0],ind, mic=True)).any():
        continue
    ind_selected += ind.tolist() # Append selected indices to list

# Positions and symbols of selected indices
pos_selected = supercell.positions[ind_selected]
symbols_selected = origsymbols[ind_selected]

# Build Atoms object excluding atoms in boundaries
finalsupercell = Atoms(
        positions=pos_selected,
        cell=supercell.cell,
        symbols=symbols_selected,
        pbc=supercell.pbc)

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