Translation Group
In this section, the main usage of the classes in translation group module will be demonstrated with examples. The module is developed in an object oriented manner, such that the Cluster and the Lattice classes serve as the foundation of the module and the other classes (e.g. Kpoints, Crystal, CrystalFTG and QpointsN) are derived from these two classes.
Cluster Class
We begin with the Cluster class, which is a collection of atoms that forms a system such as a molecule. Here we can create a Cluster class object of a \(\textrm{CO}_2\) molecule with the molecular structure illustrated in Fig. 2, and the following statements are used to construct the Cluster object for the molecule.
>>> atoms = {
"C": np.array([
[0.0, 0.0, 0.0]
]),
"O": np.array([
[ 1.0, 0.0, 0.0],
[-1.0, 0.0, 0.0]
])
}
>>> cluster = Cluster(atoms=atoms, orbitals="p_x,p_y")
Fig. 2 Structure diagram of a \(\textrm{CO}_2\) molecule.
Since the 3 atoms of the molecule are positioned along a line, the displacements can be categorized into an along-the-axis mode and a perpendicular-to-axis mode, Thus we can set the orbitals of the Cluster to \(p_x\) and \(p_y\) to simplify the study of the atomic displacements.
This class serves as a foundation to analyzing the properties of molecules, and we can apply various operations to the module, for example a \(90^\circ\) rotation:
>>> rotation_90 = np.array([
[0, 1, 0],
[1, 0, 0],
[0, 0, 1],
], dtype=float)
>>> example_cluster = cluster.copy()
>>> example_cluster.rotate(rotation_90)
>>> example_cluster.atoms
OrderedDict([('C', array([[0., 0., 0.]])), ('O', array([[ 0., 1., 0.],
[ 0., -1., 0.]]))])
Here, the center of the molecule is defined as the average positions of all atoms. In the case of the above molecule, it’s at the origin of the coordinate system:
>>> cluster.center
array([0., 0., 0.])
And we can also move the entire cluster to relocate the origin onto one of the O atoms:
>>> example_cluster = cluster.copy()
>>> example_cluster.shift_atoms([1, 0, 0])
>>> example_cluster.center
array([1., 0., 0.])
>>> example_cluster.atoms
OrderedDict([('C', array([[1., 0., 0.]])), ('O', array([[2., 0., 0.],
[0., 0., 0.]]))])
Lattice
Lattice is an infinite array of points generated by translation symmetry from the primitive translation vectors (a.k.a. lattice vectors), serving as a basic building block of the translation group.
Here the lattice vectors of face center cubic is used as a demonstration:
>>> lattice_vectors = np.array([
[0.0, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.5, 0.5, 0.0],
])
>>> lattice = Lattice(vec=lattice_vectors)
With the above Lattice class object, we can compute various properties of the lattice.
>>> print("volume:", lattice.vol)
volume: 0.25
>>> print("lengths of the lattice vectors:", lattice.abc)
lengths of the lattice vectors: [0.70710678 0.70710678 0.70710678]
>>> print("angles between the lattice vectors:", lattice.abg)
angles between the lattice vectors: [60. 60. 60.]
Meanwhile, a real space lattice comes with a reciprocal lattice. We can also compute the properties of the reciprocal lattice. (In the convention of our software, the reciprocal lattice vectors contain the \(2\pi\) factor in normalization.)
>>> print("volume of reciprocal lattice:", lattice.rvol)
volume of reciprocal lattice: 992.200854
>>> print("lengths of the reciprocal lattice vectors:", lattice.rabc)
lengths of the reciprocal lattice vectors: [10.88279619 10.88279619 10.88279619]
>>> print("angles between the reciprocal lattice vectors:", lattice.rabg)
angles between the reciprocal lattice vectors: [109.47122063 109.47122063 109.47122063]
Additionally, the Lattice class allows us to perform various operation to a lattice, for example applying an axial strain:
>>> strain = [0.02, 0, 0] # axial strain along x
>>> example_lattice = lattice.copy()
>>> example_lattice.axial_strain(strain)
>>> print(example_lattice.vec)
[[0. 0.5 0.5 ]
[0.51 0. 0.5 ]
[0.51 0.5 0. ]]
Crystal and Finite Translation Group
A crystal consists of a lattice and a collection of basis atoms in the primitive unit cell. Thus we created the Crystal class by deriving from both the Lattice and the Cluster classes simultaneously. Here we construct a rock salt crystal with lattice parameter \(a=1\) for demonstration.
>>> lattice_vectors = np.array([
[0.0, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.5, 0.5, 0.0],
])
>>> atoms = {
"Na": np.array([[0.0, 0.0, 0.0]]),
"Cl": np.array([[0.5, 0.5, 0.5]]),
}
>>> crystal = Crystal(vec=lattice_vectors, atoms=atoms)
The Crystal class inherits all the aforementioned functionalities of Lattice and Cluster classes, although the method to rotate atoms is updated due to the Crystal using direct coordinates instead of the Cartesian coordinates in Cluster class.
>>> rotation_90 = np.array([
[0, -1, 0],
[1, 0, 0],
[0, 0, 1],
], dtype=float)
>> example_crystal = crystal.copy()
>> print(example_crystal.positions_cartesian)
[[0. 0. 0. ]
[0.5 0.5 0.5]]
>> print(example_crystal.rotate_atoms(rotation_90))
[[ 0. 0. 0. ]
[ 1.5 -0.5 -0.5]]
Due to limitations of the complexity and computing power, we are unable to work with the infinite lattice. Thus we often turn to a finite translation group (see paper for details) as an approximation. A finite translation group is defined by a supercell with which there are a finite number of translation vectors that translate the primitive cell of the crystal across the space. Thus we can derive the CrystalFTG class from the Crystal class, and construct a CrystalFTG class object from a Crystal class object with a supercell matrix.
>>> supa = np.identity(crystal.dim, dtype=int) * 2
>>> supa
array([[2, 0, 0],
[0, 2, 0],
[0, 0, 2]])
>>> crystal_ftg = CrystalFTG.from_primitive(crystal, supa)
>>> crystal_ftg.to_dict()
OrderedDict([('vec',
array([[0. , 0.5, 0.5],
[0.5, 0. , 0.5],
[0.5, 0.5, 0. ]])),
('atoms',
OrderedDict([('Na', array([[0., 0., 0.]])),
('Cl', array([[0.5, 0.5, 0.5]]))])),
('orbitals', None),
('supa',
array([[2, 0, 0],
[0, 2, 0],
[0, 0, 2]]))])
The finite translation group, provides various properties of the crystal.
For example, we can construct the Wigner-Seitz cell of the supercell. This is a crucial piece of information used in our Fourier interpolation scheme (see paper for details).
>>> ws_cell_translation, ws_cell_atoms, ws_cell_weights = \
crystal_ftg.get_wigner_seitz_cell(center=[0, 0, 0])
>>> ws_cell_weights
array([1. , 0.5 , 0.5 , 0.5 , 0.25 ,
0.5 , 0.5 , 0.5 , 0.5 , 0.5 ,
0.5 , 0.5 , 0.5 , 0.5 , 0.16666667,
0.16666667, 0.16666667, 0.16666667, 0.16666667, 0.16666667,
0.25 , 0.25 , 0.25 , 1. , 1. ,
1. , 1. , 1. , 1. , 0.25 ,
0.25 , 0.25 , 0.25 ])
Additionally, we can construct the naive displacements basis vectors of a \(\textbf{q}\)-point. This is used to construct the symmetrized basis in our LID approach, where we displace atoms along either the real or imaginary part of the reciprocal space displacements basis, in order to only activate specific irreducible derivatives.
>>> qpoint = parse_array("1/2 0 0", dtype=Fraction)
>>> qpoint
[Fraction(1, 2), Fraction(0, 1), Fraction(0, 1)]
>>> crystal_ftg.orbitals = "p"
>>> basis = crystal_ftg.get_basis_at_q(qpoint)
>>> basis[(0, "p_x")].real
array([[ 0.35355339, 0. , 0. ],
[-0.35355339, 0. , 0. ],
[ 0.35355339, 0. , 0. ],
[-0.35355339, 0. , 0. ],
[ 0.35355339, 0. , 0. ],
[-0.35355339, 0. , 0. ],
[ 0.35355339, 0. , 0. ],
[-0.35355339, 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ],
[ 0. , 0. , 0. ]])
Algorithm to find translation vectors of any given supercell
Although diagonal matrices are most commonly used as supercell matrices, non-diagonal supercells are crucial in the lattice dynamic analysis of crystals. Phonons and their interactions can be extracted in the smallest non-diagonal supercells that accommodate the irreducible derivatives (see paper for details). Therefore, a robust algorithm is required to find the translation vectors of a given supercell matrix, as well as indexing the vectors with the supercell.
The details of the algorithm to find translation vectors of any given supercell is explained in the Appendix of the paper. The same algorithm can be reversed to compute the indices of the lattice points at \(\mathcal{O}(1)\) time complexity without searching through the list of all translation vectors.
>>> non_diagonal_supa = np.ones((3, 3), dtype=int) - 2 * np.identity(3, dtype=int)
>>> lattice_ftg = LatticeFTG(np.identity(3), non_diagonal_supa)
>>> lattice_ftg.lattice_points
array([[0, 0, 0],
[0, 0, 1],
[0, 1, 0],
[1, 0, 0]])
>>> lattice_ftg.get_index(lattice_ftg.lattice_points)
array([0, 1, 2, 3])
Algorithm to find the minimum supercell for a given list of \(\textbf{q}\)-points
The minimum supercell multiplicity is derived in the paper. Here using the implemented function, we find that for the \(\textbf{q}\)-point \(\textbf{q}=\left(\frac{1}{4}, \frac{3}{4}, 0\right)\) the supercell of the minimum multiplicity is:
, and for the \(Q\)-point \(Q=\left[\left(\frac{1}{4}, \frac{3}{4}, 0\right), \left(\frac{1}{4}, \frac{1}{2}, 0\right), \left(\frac{3}{4}, \frac{1}{4}, \frac{1}{2}\right)\right]\) the supercell of the minimum multiplicity is:
>>> Qpoint = parse_array("""\
1/4 3/4 0
""", dtype=Fraction)
>>> Qpoint
[Fraction(1, 4), Fraction(3, 4), Fraction(0, 1)]
>>> get_minimum_supercell(Qpoint)
array([[1, 1, 0],
[0, 4, 0],
[0, 0, 1]])
>>> Qpoint = parse_array("""\
1/4 3/4 0
1/4 1/2 0
3/4 1/4 1/2
""", dtype=Fraction)
>>> Qpoint
[[Fraction(1, 4), Fraction(3, 4), Fraction(0, 1)],
[Fraction(1, 4), Fraction(1, 2), Fraction(0, 1)],
[Fraction(3, 4), Fraction(1, 4), Fraction(1, 2)]]
>>> get_minimum_supercell(Qpoint)
array([[4, 0, 0],
[0, 4, 0],
[0, 0, 2]])
\(Q\)-points at a given order and FTG
Phonon interactions are associated with \(Q\)-points at the order of the interaction. We implemented an efficient algorithm to find symmetrically irreducible Q-points in a given FTG and order, that considers both the point symmetry and the permutation symmetry of the interaction.
Here we continue with the NaCl example. The first step is to find the irreducible \(\textbf{q}\)-points.
>>> lattice_vectors = np.array([
[0.0, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.5, 0.5, 0.0],
])
>>> supa = np.identity(crystal.dim, dtype=int) * 4
>>> pg = "Oh"
>>> qpoints = Kpoints(vec=lattice_vectors, supa=supa, pg=pg)
>>> qpoints.find_irreducible_lattice_points()
>>> qpoints.irreducible_kpoints
array([[Fraction(0, 1), Fraction(0, 1), Fraction(0, 1)],
[Fraction(1, 4), Fraction(0, 1), Fraction(0, 1)],
[Fraction(1, 2), Fraction(0, 1), Fraction(0, 1)],
[Fraction(1, 4), Fraction(1, 4), Fraction(0, 1)],
[Fraction(1, 2), Fraction(1, 4), Fraction(0, 1)],
[Fraction(3, 4), Fraction(1, 4), Fraction(0, 1)],
[Fraction(1, 2), Fraction(1, 2), Fraction(0, 1)],
[Fraction(3, 4), Fraction(1, 2), Fraction(1, 4)]], dtype=object)
At the same time we can construct the list of point group operations that transform a k-point in FTG to its irreducible counterpart, and the indices of irreducible k-point each k-point can be rotated into.
>>> qpoints.irreducible_kpoints_trans
array(['E', 'E', 'E', 'I', 'Ic2b', 'E', 'E', 'E', 'Ic2b', 'Ic2b', 'E',
'Ic2z', 'Ic2z', 'I', 'I', 'I', 'Ic2d', 'Ic2f', 'Ic2f', 'Ic2f',
'Ic2d', 'Ic2x', 'Ic3ga', 'Ic3ga', 'Ici4y', 'Ic2e', 'Ic2x', 'E',
'Ic3be', 'Ici3al', 'Ic2b', 'Ic4z', 'Ic2d', 'Ic4x', 'Ic2e',
'Ici3de', 'Ic2d', 'Ic2y', 'Ici3be', 'Ic2f', 'Ic2c', 'Ic3al',
'Ic2e', 'Ic4y', 'Ic3be', 'Ic2e', 'Ic2e', 'Ic2x', 'Ic2y', 'Ic4x',
'Ici3ga', 'Ic2e', 'Ic2d', 'Ic2y', 'Ic2a', 'Ici4x', 'Ic3de', 'I',
'Ici4z', 'Ici3be', 'Ic2c', 'Ic2c', 'Ic2c', 'Ic2e'], dtype='<U7')
>>> qpoints.irreducible_kpoints_map
array([0, 1, 2, 1, 1, 3, 4, 5, 2, 4, 6, 4, 1, 5, 4, 3, 1, 3, 4, 5, 3, 1,
5, 4, 4, 5, 4, 7, 5, 4, 7, 4, 2, 4, 6, 4, 4, 5, 4, 7, 6, 4, 2, 4,
4, 7, 4, 5, 1, 5, 4, 3, 5, 4, 7, 4, 4, 7, 4, 5, 3, 4, 5, 1])
Next, we can use the above information to find irreducible \(Q\)-points.
>>> qpointsn = QpointsN(vec=lattice_vectors, supa=supa, order=3, pg=pg)
>>> qpointsn.find_irreducible_Qpoints()
array([[[Fraction(0, 1), Fraction(0, 1), Fraction(0, 1)],
[Fraction(0, 1), Fraction(0, 1), Fraction(0, 1)],
[Fraction(0, 1), Fraction(0, 1), Fraction(0, 1)]],
...,
[[Fraction(3, 4), Fraction(1, 4), Fraction(1, 2)],
[Fraction(3, 4), Fraction(1, 4), Fraction(1, 2)],
[Fraction(1, 2), Fraction(1, 2), Fraction(0, 1)]]], dtype=object)