VASP symmetry in piezo-electric tensor

[maartendft]

Dear admin,

For compounds with (-6m2) point groups, textbooks report the symmetry of the piezo-electric tensor as:

0 0 0 0 0 -2*d_22
-d_22 d_22 0 0 0 0
0 0 0 0 0 0

When I run compounds with this symmetry through VASP, I get answers such as (in C/m^2):

0 0 0 0 0 -0.4028
-0.4028 0.4028 0 0 0 0
0 0 0 0 0 0

Note that the component in the upper right corner is missing the factor "2", compared to the symmetries reported in the literature.

My suspicion is that is has to do with the VASP-definition of the off-diagonal terms in the 3x3 block on the right. These always seem scaled by a factor 0.5 w.r.t. some other definitions. I know this sometimes also occurs when people use different definitions for the elastic tensor and strain, absorbing the factors of 2 for shear into the tensor or the strain vector.

Would it be possible to comment on this?

Thanks,

Maarten de Jong

[hipertrofia]

I would argue that the result you are getting from VASP is the correct one. In a two-dimensional hexagonal lattice with two atoms in the unit cell, for instance BN, if B is at the unit cell origin (0,0) and N is in the middle of the triangle (a0/2,sqrt(3) a0/6) then the internal strain is given by t = (xi eps6, xi [eps1-eps2]), where xi is some constant with units of Angstrom. You can check this result in Eq. (9) of PRB 86, 014117 and repercusions on the PZ tensor in Eq. (14) of PRB 88, 214103 for wurtzite.

In such a lattice the dipole moment of the unit cell depends only on internal strain (see PRB 91 075203 for the argumentation for zinc-blende), therefore dp1/deps6 = dp2/deps1 = - dp2/deps2 = xi. The polarization would then include volume terms, which happen to be trivial in the linear limit for this case (you can either use derivatives or finite differences truncated to first order terms in the strains to get this result), yielding e16 = e21 = - e22 = Z xi /V0, where V0 is the unit cell's volume and Z is some effective piezoelectric charge (in general different from the Born effective charge). For an ionic crystal, Z is the oxidation number times the elementary charge, with positive or negative sign depending on your convention for the unit cell (related through a 180 degree in-plane rotation).

In the argumentation above I have removed the third lattice vector for simplicity but with no loss of generality (it's not involved). I have a simple Mathematica script where you can obtain the internal strain result for a BN lattice. You can contact me directly if you're interested [miguel.caro (at) aalto.fi]. I'm interested in the discussion on PZ tensors anyway, so I'll be glad to receive your email.

As a final note, it is always in my opinion best to work with Voigt indices because if one does not consider rotations, then the Cartesian shear strain components are not independent and that messes up the power expansion of quantities such as the electric polarization. The result is that because of conventions and confusions these factors of 2 or 1/2 are unfortunately all over the place in the literature.

Miguel

[maartendft]

Dear Miguel,

Thanks a lot for your reply and the references, these seem very useful!

I will get back to you by e-mail once I have found time to read the papers you suggested. But before then, I'd like to mention that I did not study BN (which has a 6/mmm point group), but instead I looked at systems with the (-6m2) point group. You can check Table 2 (page 15) in the following document:

http://www.gdp.if.pwr.wroc.pl/pliki/pie ... effect.pdf

As you can see, there are different symmetry-restrictions in the piezo-electric tensor based on the 2 distinct point groups.

What VASP calculates for (-6m2) materials almost matches the tensor-symmetry as reported in the document. However, the factor "2" is missing. I suspect this is due to the factors of 2 introduced within the Voigt-notation on page 8 and 9 of the document.

Anyway, I will read your papers and get back to you about this.

Thanks,

Maarten

[hipertrofia]

The PZ tensor of bulk BN is different from what I wrote before, which is only valid for two-dimensional BN. In bulk BN the alternating hexagonal planes are rotated by 180 degrees with respect to each other, making the average PZ response zero (the response in plane A is minus the response in plane B). However this 2D case bears resemblance to your problem and so I thought it would be worth mentioning.

The document you link is very similar to the description from Nye's book (which is in the list of references) and I would assume taken from there. In Nye's book the stress components appear without the factor of 2 whereas in the linked document they appear with a factor of 2. Then, e16 = 2 e21 as given by Nye, but according to Nye's definition d16 = 2 d112. That is, there is a factor of 2 involved in going from Cartesian to Voigt PZ moduli. Again, you can grasp how problematic these convention issues are.

To be more precise, if one defines the 2D strain matrix as

eps11 eps6/2
eps6/2 eps22

and the PZ coefficients through

Pi = Sum eij epsj

then for the discussed 2D hexagonal lattice (and I'm assuming this is also true for your case) it should hold that e16 = e21 = -e22.