# Finding possible lattice planes of a crystal structure

After generating a crystal structure from a crystallographic data and duplicating it to a larger crystal system I would like to find possible lattice planes of this crystal.

It is well-known that the hexagonal ice has three major lattice planes:

• $0001$, basal face
• $10 \bar10$, 1st prismatic face
• $11\bar20$, 2nd prismatic face

but there should be more lattice planes.

So the question is how can I write a code in order to generate planes consisting of at least three local atom points.

The main idea is to construct a list of vectors between atoms (connecting atoms) and then calculate the cross product of vector pairs (two connection vectors) which are then parallel to the normal vectors of the corresponding planes.

Then a list of plane vectors should be sort by the number of equivalent plane vectors (e.g. either pointing in positive or negative x-direction, (100) or (-100)) and then the corresponding 'effective' Miller indices should be calculated.

If the crystal system of ice contains a few hundreds oxygen atoms, constructing and calculating lists of vectors would be very time-consuming. Has anyone a brilliant idea to solve this problem?

For more detail of my question above:

Here is the URL of the data file: http://www1.lsbu.ac.uk/water/ice1hsc.html (downloading the pdb file containing 1296 atoms of hexagonal ice).

After constructing a new data file containing only atom types and the (x,y,z) coordinates ot all atoms, I extract the coordinates of all oxygen atoms and make a list of these coordinates:

iceIhdata = Import["water ice Ih.dat"]
iceIhOList = Table[iceIhdata[[3*(i + 1)]], {i, 0, Length[iceIhdata]/3 - 1}].DMx[Pi/2].DMz[Pi/2]


DMx and DMz are the rotation matrices in order to get desired coordinate axes for this ice crystal. There are 432 oxygen atoms in total.

To list all oxygen-oxygen vectors:

OOvectors = DeleteDuplicates[DeleteCases[Flatten[Table[
If[i < j, iceIhOList[[i]] - iceIhOList[[j]], 0], {i, 1,
Length[iceIhOList]}, {j, 1, Length[iceIhOList]}], 1], 0]]


Then I want to calculate the normal vectors of the corresponding planes spanned by two vectors connecting two oxygens (OOvectors[i] and OOvectors[j]):

nvectors = Parallelize[Table[If[i < j, Cross[OOvectors[[i]], OOvectors[[j]]], 0], {i, 1,
Length[OOvectors]}, {j, 1, Length[OOvectors]}]]


But the PC is overheated... How can I write the code that list all vectors and count the number of equivalent planes effectively and faster?

• Ugh, that formatting was abysmal. Turns out, if you format a question properly, it suddenly becomes much clearer. Any plane is defined by three points. So you take all possible 3-tuples (each element is a point). In a lattice the atoms are regularly spaced, so if $\vec v$ connects two equivalent oxygens, $2 \vec v$ will also connect two oxygens and be part of the same plane. So you want co-prime vector components in terms of the elementary translational symmetry vectors. That's a good way to filter out very many unneeded tuples. Jun 22, 2017 at 6:21
• Voting to leave open, since it's a fun and salvageable question. Too bad, I don't have too much time on my hands today. Jun 22, 2017 at 6:22
• @LLIAMnYP: thank you for the suggestions. I will try to rewrite. Jun 22, 2017 at 7:38