[ase-users] Replacement of a tetrahedron

[Roger Mason]

Hello,

I would like to find the SiO4 tetrahedra in a silicate so that I can
replace, in one of them,  the silicon atom by a vacancy and its 4
coordinating oxygens, with OH molecules.  Ideally this should be done so
as to preserve the symmetry of the tetrahedron, so the H atoms should
lie on the vectors from the vacancy centre to the oxygen atoms.

I have been playing around with neigbor_list and borrowed heavily
(i.e. stole) from the natural_cutoff() example:

++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Origin shift +2/3 along z
q = crystal(('Si','O'), basis=[(0.46970, 0.00000, 2/3),(0.41350,0.26690,0.785766666667)], spacegroup=154, cellpar=[4.916,4.916,5.4054,90,90,120])# * (2,2,2)

# Adapted from example here:
# https://wiki.fysik.dtu.dk/ase/ase/neighborlist.html?highlight=neighbors#ase.neighborlist.natural_cutoffs
#cutOff = neighborlist.natural_cutoffs(q)
cutOff = [2.77, 2.77, 2.77, 0.66, 0.66, 0.66, 0.66, 0.66, 0.66]
#print(cutOff)
neighborList = neighborlist.NeighborList(cutOff, self_interaction=False, bothways=False)
neighborList.update(q)
matrix = neighborList.get_connectivity_matrix()
n_components, component_list = sparse.csgraph.connected_components(matrix)
idx = 1
molIdx = component_list[idx]
print("There are {} molecules in the system".format(n_components))
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

The result:

There are 1 molecules in the system

Clearly the individual tetrahedra are not being found.

Can someone please help me out?

Thanks,
Roger