big error in C16 of HCP/wurtzite strucutre

Dear lammps users and developers,
    I calculated elastic constants of wurtzite CdSe using elastic script in the examples dir. C16 I got is around 4.05 GPa, which is very strange since its principle value should be 0. this value also converged very well at different deformation sizes, so it is not a convergence error, i think.
    Did you ever calculate elastic constants of HCP or wurtzite structures using this script? I don't know if such a big error also appear in other wurtzite or HCP materials. The scripts in my calculation are also attached, any suggestions would by highly appreciated! Thanks!

Best regards,

data.lmp (505 Bytes)

displace.mod (3.38 KB)

in.elastic (5.62 KB)

init.mod (752 Bytes)

potential.mod (712 Bytes)

I believe people have characterized elastic constants for HCP structures. Please note that an empirical potential you use does not guarantee good properties compared with electronic structure methods and/or experimental values. It is just an empirical potential fitted to some target. You first have to look at the publication to see if C16 was even included in the fitting target, and if it did, see how well it reproduces C16.


Ning may be running into a different issue here. The C16 matrix element for HCP should be 0 due to symmetry reasons. The value reported by the script results from symmetrization of the 16 and 61 elements of the matrix which turn out to be quite different when printed independently, i. e., one about 8GPa and the other ~0GPa. This does not happen to the 12 and 21 elements which individually seem to preserve the symmetric character. Aidan may be the best person to ask as I believe he wrote that script.

I confirmed what Carlos found. Probably this is due to the structure
*not* being wurtzite, or even a crystal. The initial relaxed value of
pxy is non-zero, which is a red flag. I repalced the data file with a
wurtzite structure generated by LAMMPS and got the following elastic
constants first time:

Elastic Constant C11all = 65.9336257924683 GPa
Elastic Constant C22all = 65.9688151923895 GPa
Elastic Constant C33all = 64.3984797548074 GPa
Elastic Constant C12all = 37.9615565299385 GPa
Elastic Constant C13all = 32.8869534343692 GPa
Elastic Constant C23all = 32.9273899804094 GPa
Elastic Constant C44all = 15.2420561304232 GPa
Elastic Constant C55all = 15.2420561774616 GPa
Elastic Constant C66all = 13.9558365418157 GPa
Elastic Constant C14all = -1.31356050994215e-08 GPa
Elastic Constant C15all = -1.04660388610634e-08 GPa
Elastic Constant C16all = 7.54866371667155e-08 GPa
Elastic Constant C24all = -3.22779657858733e-08 GPa
Elastic Constant C25all = 1.0476101624966e-08 GPa
Elastic Constant C26all = -8.62402675904665e-08 GPa
Elastic Constant C34all = -5.5065623543402e-09 GPa
Elastic Constant C35all = 7.54523129362618e-09 GPa
Elastic Constant C36all = -3.54005816422945e-08 GPa
Elastic Constant C45all = -1.37152968761703e-08 GPa
Elastic Constant C46all = 2.04509223679242e-08 GPa
Elastic Constant C56all = 3.40997023815799e-08 GPa

the init.mod file is:

# NOTE: This script can be modified for different atomic structures,
# units, etc. See in.elastic for more info.

The Q felt legit enough that I passed it on without checking the initial struct. Only had a unix terminal anyways and little time thus skipped visualization completely. Good to know the score remains: lammps 1 user 0

Dear Aidan, Carlos and Ray,
   Thanks very much for the explanation of Carlos and confirmation of Aidan. After carefully checking my relaxed structure and that of Aidan, I found we didn't catch the key point and it is not due to the initial structure. Aidan used an 8-atom orthogonal unit cell with the perfect wurtzite lattice as his initial structure, and in my calculation is a 4-atom typical primitive cell for wurtzite. In Aidan's elastic script, structures will be relaxed first and then deformed to generate stresses. Actually, the two relaxed structures are exactly the same. please visualize them carefully.
   I guess the problem may be due to the shape of the unit cell, but cannot offer a deeper insight. Anyway, C11, C12, C13 C33 and C44 are reliable no matter which kind of unit cell is used, which is really encouraging.
   Thanks again!

Best regards,

The Cij computation should be independent of the unit cell chosen. Its been a while since I dealt with this topic but if I recall correctly quasi-static deformations (strains) for elastic constant calculations have to be symmetrized to get Cijs that do not depend of the cell shape. I cannot really comment more because I have not spent time understanding the key points of Aidan’s script.

Ning Wang,

You are right. There is something wrong with how C16 is calculated in
non-orthorhombic cells. The problem shows up for any hexagonal cell,
but not for the equivalent orthorhombic cell. When the cell is
strained in y (normal to the 6-fold axis), the pxy component relaxes
to zero under internal atomic displacements, making C62 zero.
however, when the cell is strained axially in x (along the 6-fold
axis) the pxy component of the stress tensor does not relax to zero,
making C16 non-zero. C61, C62, C16, and C26 are averaged in the
script to give a non-zero value for C16all. Except for C16, no other
elastic constants seem to be affected.


I have made some inquiries with crystal structure experts, so far no
definitive response. Purely empirically, I have observed that the
calculated C16 can be flipped from 4.05 GPa to -4.05 GPa by flipping
the sign of the xy 'tilt factor' in data.lmp. This corresponds to
flipping the cell from a forward-leaning parallelogram to an
equivalent backwards-leaning one, and so all other properties of the
crystal are unaffected. Taking the arithmetic average of these two
cases yields the desired result C16=0. It would be nice to be able to
generalize this in a more formal way and incorporate it into the
script or at least document it.

Been trying to look into it but no much time on this side. Far from being a Xtal struct expert yet my physical intuition cannot cope with the idea of lattice dependent Cijs…
Your empirical observation related to the ± changes feels as if the issue somehow has to do with symmetrization… could it be related to the fact that Lammps uses upper triangular algebra to represent the matrix tensor? As a complementary test I would try to find another code (GULP comes to mind) to run the same experiment.

one point that comes to my mind is the following:

most codes where you can input primitive cells (e.g. plane wave pseudopotential DFT codes), use different orientations in space than LAMMPS, to increase symmetries. LAMMPS however requires axes to be aligned in a particular way. that would also have consequences as to how a deformation is x-, y- or z- direction is performance relative to the orientation of the cell.


Thats what I meant by upper triangular algebra for the h matrix.

I found the problem. It was even simpler than the upper triangular
issue, I am embarassed to say. When applying axial strain in x, my
script was changing lx, but not xy, so the strain was not pure axial.
The modification to the script is easy. Now I am testing it on the
diamond primitive cell, which is a little more complicated to
construct than I expected. I will post the corrected script once it
passes that test.

I have checked in a modified script that eliminates the bad value for
C16 in Ning Wang's example.

Thanks Aidan for all your work and thanks to Ning for circumventing our earlier criticisms while keeping insisting on the issue.

Hi Aidan

I am glad that you found the problem with your script, but I just wanted to comment on something that you mentioned earlier in you thread about flipping the sign of the tilt-factor changing the sign of C16. While that was due to a bug in your script, the issue that you hinted at is real and indeed an issue for some systems. The sign of the off-diagonal elastic constants in hexagonal or rhombohedral (and perhaps other symmetries) depends on the setting of the cell.

We came across a classic example of this some years ago where computation could not reproduce the experimental elastic constants of alumina (Al2O3) — though the magnitudes were similar the sign differed. Through a careful experiment, where we took great care to track the orientation of the crystal, we reconciled theory and experiment, and in the process showed convincingly that the original, definitive experiment at NBS (now NIST) was incorrect due to presumably turning the crystal around. You can find the paper detailing this at

Almost all experimental determinations of the elastic constants propagate this error, presumably because experiment easily gives the magnitude but not the sign, so people simply followed the sign of the NBS work.

I just wanted to alert you and the readers of the list to this issue. It is perhaps not of great practical importance but it is good to be aware that such issues do exist, and that unless you are very careful, the structure that you build in LAMMPS may correspond to a non-standard setting of the crystal.

Hi Paul and Aidan,
If you allow me I will attempt to write a quick summary for future readers because the threat has turned long
and people could at this point easily get lost or take some past comments out of context.
I must add that I feel partially responsibility for doing this given that it was my original analysis on Ning’s
anomaly observation the one that triggered Aidan’s review.
To all those people who will look into this threat I think it is safe to say that the sign of certain elements
of the elastic constant tensor can depend on the orientation of the lattice. However, the properties of the tensor,
i.e., absolute magnitudes of the Cijs and their symmetric character cannot depend on which unit cell is chosen
to represent the material. Different cell geometries corresponding to the same solid should generate the same
elastic tensor (up to a sign difference). Lack of symmetry in the tensor matrix (Cij != Cji) should be seen as a sign
of trouble, and could be due to non-properly-applied lattice strains and/or incomplete atomic relaxations that leave
residual stresses acting on the cell.
I believe this covers it all but you or any other user/developer is welcome to add more info at any time.

Below always meant to say thread no threat. The price for composing messages at the exhaustion threshold.


I agree with your synopsis.