A Script for angle optimization

# 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.

SQS; elasticity; projection; tensor; rotation

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

 Ccub≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝C11C12C12000C12C11C12000C12C12C11000000C44000000C44000000C44⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (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 alTasná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):

 Ctric≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝43616116012112545316041514281338188129\omit\span\omit\span\omitSYM1869189⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (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

 Ccub=PcubCtric, (3)

which yields average projected cubic elastic constants of

 ¯Ccub11=Ctric11+C% tric22+Ctric333=439.0GPa, ¯Ccub12=Ctric12+C% tric13+Ctric233=160.3GPa, ¯Ccub44=Ctric44+C% tric55+Ctric663=187.7GPa, (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)

 Csym=N∑i=1aiVi, (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

 R(θx,θy,θz)≡Rz(θz)Ry(θy)Rx(θx). (6)

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

 Rx=⎛⎜⎝1000cosθx−sinθx0sinθxcosθx⎞⎟⎠, Ry=⎛⎜⎝cosθy0sinθy010−sinθy0cosθy⎞⎟⎠, Rz=⎛⎜⎝cosθz−sinθz0sinθzcosθz0001⎞⎟⎠. (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

 ~C(θx,θy,θz)≡R(θx,θy,θz)C. (8)

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

 ~Cijkl=∑m,n,o,pRimRjnRkoRlpCmnop, (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

 Csym(θx,θy,θz)=Psym~C(θx,θy,θz). (10)

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

 ∂∂ai∥~C(θx,θy,θz)−Csym(θx,θy,θz)∥2=0, (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:

 ∂∂θi∥~C(θx,θy,θz)−Csym(θx,θy,θz)∥2=0. (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

 ∥~C(θx,θy,θz)−P% sym~C(θx,θy,θz)∥2 (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 amount in the range between and  GPa to each component of the perfect hexagonal elastic tensor of wurtzite GaN (in GPa): Caro et al. (2012)

 C≡⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝352155852−19−83789413−8−25395−24−18−1910305\omit\span\omit\span\omitSYM11115118⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (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:

 θcubx≈ C34−C24Dcub,θcuby≈C15−C35Dcub, θcubz≈ C26−C16Dcub, Ccub11≈ 13(C11+C22+C33)+43E% cub, Ccub12≈ 13(C12+C13+C23)−23E% cub, Ccub44≈ 13(C44+C55+C66)−23E% cub, (15)

with

 Dcub= 23(C11+C22+C33)−23(C12+C13+C23) −43(C44+C55+C66), Ecub= (C34−C24)θcubx+(C15−C35)θcuby +(C26−C16)θcubz. (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

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

### References

1. S.-H. Wei, L. G. Ferreira, J. E. Bernard,  and A. Zunger, Phys. Rev. B 42, 9622 (1990).
2. P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
3. W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
4. F. Tasnádi, M. Odén,  and I. A. Abrikosov, Phys. Rev. B 85, 144112 (2012).
5. M. Moakher and A. N. Norris, J. Elasticity 85, 215 (2006).
6. Ç. Diner, M. Kochetov,  and M. A. Slawinski, J. Elasticity 102, 175 (2011).
7. True random numbers have been obtained from www.random.org.
8. M. A. Caro, S. Schulz,  and E. P. O’Reilly, Phys. Rev. B 86, 014117 (2012).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minumum 40 characters