# Extended scheme for the projection of material tensors of arbitrary symmetry onto a higher symmetry tensor

## Abstract

I propose a straightforward generalization of the projection scheme for elastic tensors introduced by Moakher and Norris [J. Elasticity 85, 215 (2006)] that takes into account also rotations. The “closest” tensor of any desired symmetry to the original tensor of lower symmetry is “closer” in this generalized scheme. The method has an important application in the context of the special quasirandom structure (SQS) method for the computational modeling of alloys, whereby the supercell’s symmetry, and therefore that of the tensors representing its properties, is reduced with respect to the material’s underlying symmetry. The approach allows to extract the tensor components most representative of the macroscopic symmetry of the material. Although the approach is general, in the present case I apply it to the elastic tensor and give numerical examples. Simple approximate analytical expressions for cubic materials are also provided.

In the context of material science, many material properties
are represented by tensors of different ranks. The electric
polarization is given by a rank-1 tensor , strain is
represented by a rank-2 tensor , the piezoelectric
response is given by a rank-3 tensor , elasticity
can be described with a rank-4 tensor , and so on.
The geometrical particularities of the different materials give
rise to the different symmetries of the tensor representing
their properties, as determined by the material’s point
group: hexagonal, cubic, etc.
In some cases, these symmetries are “slightly” broken, and
materials with an underlying expected symmetry might present
deviations due to imperfections, impurities, alloying effects,
and similar. Within the frame of computational materials science
it is common practice, in part due to computational limits,
to use finite size supercells to represent actual materials.
The special quasirandom structure (SQS) method Wei *et al.* (1990)
is a common approach when modeling alloys in the context
of ab initio calculations, such as density functional theory
(DFT). Hohenberg and Kohn (1964); Kohn and Sham (1965) The use of
these SQS leads to supercells whose properties slightly deviate
from the ones expected from the crystallographic class of the
macroscopic alloy. These are due to the small deviations of the finite-size
cell with respect to the perfect supercell geometry of the macroscopic
material.
Take as an example the case of the stiffness
tensor of a cubic material which, employing Voigt notation
for the sake of convenience, can be written in matrix form as

(1) |

where the elements (with uppercase indices indicating Voigt notation) are the elastic constants of the material. There are two sets of symmetries applicable to . The first one is the general intrinsic symmetry of the stiffness tensor upon exchange of . The second one arises from the specific cubic symmetry of the crystal class, whereby all three Cartesian directions are equivalent within this preferential reference frame (I shall come back to this later on), and therefore , and .

When using the SQS approach, e.g. to study the elastic properties of
a cubic alloy, one is left using a “pseudo-cubic” unit cell which could have
its symmetry reduced to anything lower than cubic, including
triclinic (i.e. no symmetry at all). Tasnádi
et al. Tasnádi *et al.* (2012) used the SQS approach together
with Moakher and Norris’s projector scheme, Moakher and Norris (2006) that I shall
discuss in more detail in the following, to study cubic TiAlN
alloys. One of the triclinic supercells gave the
following stiffness tensor (in GPa):

(2) |

which resembles a cubic stiffness tensor [Eq. (1)] but strictly speaking has
inherited the triclinic symmetry (no symmetry) of the parent supercell. In
order to extract the cubic “part” of Eq. (2), Tasnádi et al.
resorted to the projection scheme introduced by Moakher and Norris. This
method relies on projector operators
that project any given tensor,
for instance in Eq. (2),
onto the closest tensor of a higher
symmetry of choice , by means of minimizing the
Euclidean distance between the two, in this case. The reader is referred
to Ref. Moakher and Norris, 2006 for the details of the method.
Moakher and Norris’s approach makes for an elegant and powerful formalism.
In the particular case of Eq. (2),
the projection onto
is done using (see
Refs. Tasnádi *et al.*, 2012; Moakher and Norris, 2006 for its specific form) as

(3) |

which yields average projected cubic elastic constants of

(4) |

and Euclidean distance .
The problem with this approach is that it does not take into account the
rotational degrees of freedom. The cubic elastic tensor takes the form
indicated in Eq. (2) only if the crystallographic axes coincide with the
Cartesian axes in which the representation of the tensor is carried out.
Because the material properties are not affected by rotations, one can
further reduce the Euclidean distance between
and by rotating either one of them so as to
maximize the projection. In fact, Moakher and Norris’s projector scheme
implemented as is would lead to the inconsistency that a perfectly
cubic elastic tensor, rotated with respect to the reference frame in which
the projector is obtained, would differ from its own
cubic projection in that
reference frame. The issue of system orientation was already highlighted in
their original work. Moakher and Norris (2006)
If the triclinic stiffness tensor of Eq. (2) is rotated
by , and degrees (where “” means
counterclockwise and
“” clockwise) with respect to the first, second and third
Cartesian axes respectively, in that order, then the projection
of onto its closest cubic elastic tensor
yields ,
and
. More importantly, the
Euclidean distance between
and is further reduced to 83.7 GPa, where
the tilde denotes the rotation performed. Note that in this case the
difference between the two different approaches is small because
was already “almost” cubic (the rotation
angles to correct the structure are correspondingly small). In a more general
case, the input tensor might be in a form that does not closely
resemble the symmetry of interest, leading to even
larger discrepancies. Based on similar considerations, Diner et al.
have recently presented a rotation-based method to identify the symmetry
class of a tensor which is not expressed in its natural coordinate
system. Diner *et al.* (2011)

Since, as previously mentioned, rotations do not affect the properties of materials, only their mathematical representation, the procedure introduced here allows to obtain a better (“closer”) higher-symmetry projection of the original tensor. In the following I deal with the details of the present approach in the context of Moakher and Norris’s method, whose details are given in Ref. Moakher and Norris, 2006.

Let and be material tensors of arbitrary symmetry and specific symmetry “sym”, respectively. Then can be expressed as a linear combination of basis tensors as Moakher and Norris (2006)

(5) |

where are constant coefficients and is the size of the basis. For the elasticity of cubic materials, there are 3 independent elastic constants and therefore .

Now, as I have discussed, a rotation of does not affect its properties, only its representation, and thus I can define the general rotation operator

(6) |

In matrix form, each individual rotation takes the following form:

(7) |

The different give counterclockwise rotation angles around the Cartesian axes. Note that the order in which the different rotations are carried out is important in determining . Therefore the general form of tensor is given by

(8) |

In the particular case of the stiffness tensor (rank-4), the matrix elements of are given by

(9) |

where all the indices run from 1 to 3. I could have equivalently made the rotation operator act on the different basis components , however this would lead to complicated angle-dependent projectors. It seems best and simplest to keep Moakher and Norris’s original formulation for the projectors and introduce the angular dependence on instead. The projector (see Ref. Moakher and Norris, 2006) for a given symmetry sym links and as

(10) |

The condition that the Euclidean distance between and be minimized is given by

(11) |

which is the original condition in the formulation of Moakher and Norris, Moakher and Norris (2006) plus a new requirement for the rotation angles:

(12) |

With the first condition only, Eq. (11), the projector can be obtained analytically from the alone, independent of , as done by Moakher and Norris. This can then be fed into Eq. (12) and the quantity

(13) |

minimized numerically with respect to the rotation angles. Alternatively, it can be shown that this minimization is equivalent to the condition of maximum for the Euclidean norm of the projected tensor, , taking the rotation angles as variational parameters. To carry out this task in the previous example, I have used the Mathematica symbolic calculator, using the same script as given in the Appendix. The procedure is straightforward and can be readily extended to other computational tools.

Note that if an isotropy plane is present then one of the rotation angles
becomes redundant. For example, a hexagonal projection which takes the
axis as parallel to maintains a constant Euclidean distance
with the original triclinic tensor for any arbitrary value of .
This allows to effectively visualize the effect of the present approach
by plotting the Euclidean distance for a hexagonal projection as a
function of and . In order to do this I have
generated the following triclinic elastic tensor by adding a
random ^{1}*et al.* (2012)

(14) |

Figure 1 shows a plot of the Euclidean distance between the triclinic elastic tensor given in Eq. (14) and its hexagonal projection as a function of and , which as discussed is independent of . It can be observed that a careless projection without considering the rotational degrees of freedom would lead to unoptimized calculated elastic constants, projected at . The optimal hexagonal projection requires a previous rotation corresponding to approximately , .

In the limit where the triclinic tensor resembles the symmetry expected, for instance cubic in the case of Eq. (2), one can make a number of assumptions that allow to obtain analytical expressions for the projected elastic constants. In particular, one can make the assumptions that i) the components of the rotated triclinic tensor are linear in the rotation angles ( and for small ), ii) that the off-diagonal components of the triclinic tensor are small compared to the block-diagonal ones (e.g. , etc.) and iii) that the difference between the components that should be equal by symmetry for the projected tensor is small, for instance , etc. in the cubic case. For the cubic crystal class, the resulting approximate expressions are:

(15) |

with

(16) |

Equivalent expressions for other crystal classes are considerably lengthier than for the cubic one, and therefore I have opted to only report the latter as way of example. Obviously, the numerical solution is always more accurate, in particular when large rotation angles are required and/or the tensor does not closely resemble the target projection symmetry.

In summary, I have presented an extension to Moakher and Norris’s formalism Moakher and Norris (2006) that allows to find the material tensor closest to a tensor of lower symmetry including also the rotational degrees of freedom. These rotational degrees of freedom do not determine the symmetry properties of the tensor but rather their mathematical representation and therefore ought to be considered when searching for higher symmetry representations of the tensor’s properties. I have explicitly carried out this calculation for the stiffness tensor and exemplified it in the context of the SQS approach, for which these considerations are readily applicable. Approximate analytical expressions have been provided for the case of the cubic crystal class.

The author is thankful to Rémi Zoubkoff for the critical reading of this manuscript.

*

## Appendix A Script for angle optimization

The following Mathematica script performs the angle-optimized projection of the triclinic elastic tensor given in Eq. (2) leading to the calculation of the projected elastic constants and rotation angles given throughout the text.

(* Declare triclinic elastic tensor and empty rotated triclinic tensor *) ctensor = Array[ct, {3, 3, 3, 3}]; rotctensor = Array[0 &, {3, 3, 3, 3}]; (* Define rotation matrices *) Rx = {{1, 0, 0}, {0, Cos[tx], -Sin[tx]}, {0, Sin[tx], Cos[tx]}}; Ry = {{Cos[ty], 0, Sin[ty]}, {0, 1, 0}, {-Sin[ty], 0, Cos[ty]}}; Rz = {{Cos[tz], -Sin[tz], 0}, {Sin[tz], Cos[tz], 0}, {0, 0, 1}}; R = Rz.Ry.Rx; (* Rotate triclinic elastic tensor *) Do[rotctensor[[i, j, k, l]] = Sum[R[[i, m]]*R[[j, n]]*R[[k, o]]*R[[l, p]]*ctensor[[m, n, o, p]], {m, 1, 3}, {n, 1, 3}, {o, 1, 3}, {p, 1, 3}], {i, 1, 3}, {j, 1, 3}, {k, 1, 3}, {l, 1, 3}]; (* Assign numerical values to the elements of the triclinic elastic tensor. Note symmetries applied *) ct[1, 1, 1, 1] = 436; ct[2, 2, 2, 2] = 453; ct[3, 3, 3, 3] = 428; ct[1, 1, 2, 2] = ct[2, 2, 1, 1] = 161; ct[1, 1, 3, 3] = ct[3, 3, 1, 1] = 160; ct[2, 2, 3, 3] = ct[3, 3, 2, 2] = 160; ct[1, 2, 1, 2] = ct[1, 2, 2, 1] = ct[2, 1, 1, 2] = ct[2, 1, 2, 1] = 189; ct[1, 3, 1, 3] = ct[1, 3, 3, 1] = ct[3, 1, 1, 3] = ct[3, 1, 3, 1] = 186; ct[2, 3, 2, 3] = ct[2, 3, 3, 2] = ct[3, 2, 2, 3] = ct[3, 2, 3, 2] = 188; ct[1, 1, 1, 2] = ct[1, 1, 2, 1] = ct[1, 2, 1, 1] = ct[2, 1, 1, 1] = 25; ct[1, 1, 1, 3] = ct[1, 1, 3, 1] = ct[1, 3, 1, 1] = ct[3, 1, 1, 1] = 11; ct[1, 1, 2, 3] = ct[1, 1, 3, 2] = ct[2, 3, 1, 1] = ct[3, 2, 1, 1] = 12; ct[2, 2, 2, 1] = ct[2, 2, 1, 2] = ct[2, 1, 2, 2] = ct[1, 2, 2, 2] = 1; ct[2, 2, 2, 3] = ct[2, 2, 3, 2] = ct[2, 3, 2, 2] = ct[3, 2, 2, 2] = 4; ct[2, 2, 1, 3] = ct[2, 2, 3, 1] = ct[1, 3, 2, 2] = ct[3, 1, 2, 2] = 15; ct[3, 3, 3, 2] = ct[3, 3, 2, 3] = ct[3, 2, 3, 3] = ct[2, 3, 3, 3] = 13; ct[3, 3, 3, 1] = ct[3, 3, 1, 3] = ct[3, 1, 3, 3] = ct[1, 3, 3, 3] = 3; ct[3, 3, 2, 1] = ct[3, 3, 1, 2] = ct[2, 1, 3, 3] = ct[1, 2, 3, 3] = 8; ct[1, 2, 1, 3] = ct[1, 2, 3, 1] = ct[2, 1, 1, 3] = ct[2, 1, 3, 1] = ct[1, 3, 1, 2] = ct[1, 3, 2, 1] = ct[3, 1, 1, 2] = ct[3, 1, 2, 1] = 9; ct[2, 1, 2, 3] = ct[2, 1, 3, 2] = ct[1, 2, 2, 3] = ct[1, 2, 3, 2] = ct[2, 3, 2, 1] = ct[2, 3, 1, 2] = ct[3, 2, 2, 1] = ct[3, 2, 1, 2] = 9; ct[3, 2, 3, 1] = ct[3, 2, 1, 3] = ct[2, 3, 3, 1] = ct[2, 3, 1, 3] = ct[3, 1, 3, 2] = ct[3, 1, 2, 3] = ct[1, 3, 3, 2] = ct[1, 3, 2, 3] = 12; (* Express the rotated triclinic elastic tensor in vector form, with 21 components *) (* Note the Sqrt[2], 2 and 2 Sqrt[2] factors to preserve the norm *) rotcvector = {rotctensor[[1, 1, 1, 1]], rotctensor[[2, 2, 2, 2]], rotctensor[[3, 3, 3, 3]], Sqrt[2] rotctensor[[2, 2, 3, 3]], Sqrt[2] rotctensor[[1, 1, 3, 3]], Sqrt[2] rotctensor[[1, 1, 2, 2]], 2 rotctensor[[2, 3, 2, 3]], 2 rotctensor[[1, 3, 1, 3]], 2 rotctensor[[1, 2, 1, 2]], 2 rotctensor[[1, 1, 2, 3]], 2 rotctensor[[2, 2, 1, 3]], 2 rotctensor[[3, 3, 1, 2]], 2 rotctensor[[3, 3, 2, 3]], 2 rotctensor[[1, 1, 1, 3]], 2 rotctensor[[2, 2, 1, 2]], 2 rotctensor[[2, 2, 2, 3]], 2 rotctensor[[3, 3, 1, 3]], 2 rotctensor[[1, 1, 1, 2]], 2 Sqrt[2] rotctensor[[1, 3, 1, 2]], 2 Sqrt[2] rotctensor[[2, 3, 1, 2]], 2 Sqrt[2] rotctensor[[2, 3, 1, 3]]}; (* Obtain the cubic projector following Moakher and Norris' recipe *) (* Define the three basis elements for the cubic elastic tensor in vector form *) velacub[1] = {1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; velacub[2] = {0, 0, 0, Sqrt[2], Sqrt[2], Sqrt[2], 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; velacub[3] = {0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; (* Generate projector *) delacub = Array[0 &, {3, 3}]; Do[delacub[[i, j]] = velacub[i].velacub[j], {i, 1, 3}, {j, 1, 3}] idelacub = Inverse[delacub]; Pelacub = Sum[idelacub[[i, j]]*Outer[Times, velacub[i], velacub[j]], {i, 1, 3}, {j, 1, 3}]; (* Numerical minimization of the Euclidean distance between the original triclinic tensor and its cubic projection with respect to the rotation angles. Optimized angles tx, ty and tz are in radians *) NMinimize[(rotcvector - Pelacub.rotcvector).(rotcvector - Pelacub.rotcvector), {tx, ty, tz}] {6999.66, {tx -> -0.0329499, ty -> -0.0319465, tz -> 0.111128}} (* Evaluate cubic projection of the rotated tensor for those angles. Gives the 21 components of the elastic tensor *) N[Pelacub.rotcvector] /. tx -> -0.0329499 /. ty -> -0.0319465 /. tz -> 0.111128 {436.836, 436.836, 436.836, 228.276, 228.276, 228.276, 377.497, 377.497, 377.497, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

When the normalizing factors are taken into account ( for and 2 for ) the elastic constants of the closest cubic projection of the original triclinic elastic tensor are obtained.

### Footnotes

- True random numbers have been obtained from www.random.org.

### References

- S.-H. Wei, L. G. Ferreira, J. E. Bernard, and A. Zunger, Phys. Rev. B 42, 9622 (1990).
- P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
- W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
- F. Tasnádi, M. Odén, and I. A. Abrikosov, Phys. Rev. B 85, 144112 (2012).
- M. Moakher and A. N. Norris, J. Elasticity 85, 215 (2006).
- Ç. Diner, M. Kochetov, and M. A. Slawinski, J. Elasticity 102, 175 (2011).
- True random numbers have been obtained from www.random.org.
- M. A. Caro, S. Schulz, and E. P. O’Reilly, Phys. Rev. B 86, 014117 (2012).