Question on elastic constants calculation using Example code

Dear LAMMPS users,

Hi,

Using Elastic constants calculation codes available in LAMMPS directory I got the attached result. C11, C12, C33, and C13 has been calculated and reported experimentally and The result I obtained are in excellent agreement with them. However, C16 C26 C36 are not obtained zero from the code. While I personally expected them to be zero. So, I wonder if you could please share some comments with me about it.

Regards,

Ali

Result.PNG

Ali,

There really isn’t enough detail here to make anything more than a guess. You have a material with cubic symmetry I am guessing? But which potential? Crystal structure? Lattice constant? Displacement magnitude? LAMMPS version? The list goes on.

-Mitch

Hi Mitch,

Thanks for the response. For the potential, I’m using a tabulated potential, which is not available in LAMMPS directory. It has a monoclinic lattice unit. Displacement is 1E-6 and the result is not sensitive to it. I wrapped the non-cubic unit cell into a cubic one. And, I’m using the latest Sep 1st 2017 version. I checked the stress component and found out that for C16 for instance, the simulation gives me a high value of Stress xy upon stretching in x direction which cause that high C16 !!!
Thanks,
Ali

As always, we suggest you use the most recent version of LAMMPS which is currently Mar 13, 2018. When I say cubic symmetry, I mean the lattice and basis that defines your material, not the shape of the simulation cell. Translating a non-cubic to a cubic simulation cell will result in some lattice strain, even if done correctly. I would suggest using the original tilted simulation cell for your tests, not the orthorhombic cell that you forced it into (liming the sources of possible error). But your material is monoclinic, why should you expect that C_*6 be zero?

“But your material is monoclinic, why should you expect that C_*6 be zero? “

Early morning mistake on my part, I was confusing mono-clinic with tri-clinic. Only C66 should be non-zero in your case.

Hi Mitch,

That’s correct. When by updating the LAMMPS I got desired result. I also have a question and I would really appreciate in case you can comment on it. I’m working on a ceramic system. Based on my simulations, I got two interesting point:

  1. The code uses C11 = 0.5*( C11-posetive + C11-negative). For metallic system Those C11-pos and C11-neg should be the same. But, Is it so for ceramic as well ? I asked because ceramics are brittle and for example for my system I got C11-neg = 1009GPa and C11-pos = 12 which means C11= 510.5 a value in good agreement with experimental data.
  2. I attached an image showing the distribution of atomic stress for my unit cell. As you can see, Si ( Blue ones ) atoms are under tension, while N ( Red ones ) atoms are under compression. But, the total pressure of the system is around 0.

What do you think of these two points?

Thanks in advance,

Ali

Capture.PNG

Ali,

  1. I would expect a little difference(few GPa if the displacement is chosen correctly) between the pos/neg results, but 1TPa differences seems wrong to me.

  2. is quite an interesting result, I assume the stress units are GPa which are quite large for an equilibrated structure! Are you sure this phase of Si-N is the predicted ground state for this interatomic potential?

-Mitch