[ase-users] Replacement of a tetrahedron

[Keilbart, Nathan Daniel]

Hello Roger,

If I remember from my crystal chemistry course, you can multiply the natural cutoffs by up to 1.15 to find all nearest neighbors that might create bonds. Additionally, I would recommend using the ConvexHull function from scipy.spatial. I then check to make sure that the central atom, Si in your case, does not lie on the hull. From there you should know which atoms are at the center of a tetrahedra and can replace as you please. I hope this has helped.

Nathan

-----------------------------------------------------------------------------------
Nathan Keilbart, PhD
Postdoctoral Research Scientist, Quantum Simulations Group
Lawrence Livermore National Laboratory
(925) 423-6620
-----------------------------------------------------------------------------------

Message: 1
Date: Sun, 26 Apr 2020 09:11:36 -0230
From: Roger Mason <rmason at mun.ca<mailto:rmason at mun.ca>>
To: ASE List <ase-users at listserv.fysik.dtu.dk<mailto:ase-users at listserv.fysik.dtu.dk>>
Subject: [ase-users] Replacement of a tetrahedron
Message-ID: <y65v9lmjwy7.fsf at mun.ca<mailto:y65v9lmjwy7.fsf at mun.ca>>
Content-Type: text/plain

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


------------------------------

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20200427/0ad0f9d4/attachment-0001.html>