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?