Elastic Constants at higher Temperature

[Habib]

Hi there,

I'm currently working on simulating the Elastic Constants of TiN, as attached, at various temperatures ranging from 0 K to 1500 K. DFT simulations allow us to calculate these constants only at 0 K. There are several approximations available to estimate the Elastic Constants at elevated temperatures https://doi.org/10.1016/j.tsf.2021.138872. I am using QuantumATK Q-2019.12. Could you please provide guidance, scripts, or options within QuantumATK that could assist me in simulating the Elastic Constants of crystal materials at different temperatures?

Elastic constants of TiN unit cell at 0K
+------------------------------------------------------------------------------+
| Elastic Constants in GPa                                                     |
+------------------------------------------------------------------------------+
|   523.42     136.72     136.72       0.00       0.00       0.00              |
|              523.42     136.72       0.00       0.00       0.00              |
|                         523.42       0.00       0.00       0.00              |
|                                    138.84       0.00       0.00              |
|                                               138.84       0.00              |
|                                                          138.84              |
+------------------------------------------------------------------------------+
+------------------------------------------------------------------------------+
| Material properties calculated from the elastic constants:                   |
+------------------------------------------------------------------------------+
| Moduli in units of GPa:                                                      |
+------------------------------------------------------------------------------+
|                 Reuss     Voigt     Hill                                     |
+------------------------------------------------------------------------------+
| Bulk modulus:   265.6166  265.6166  265.6166                                 |
| Shear modulus:  156.4882  160.6448  158.5665                                 |
+------------------------------------------------------------------------------+
|                     X         Y         Z                                    |
| Young's modulus:   466.7859  466.7859  466.7859                              |
+------------------------------------------------------------------------------+
|                     XY        XZ         YZ                                  |
| Poisson ratios:      0.2071    0.2071    0.2071                              |
|                     YX        ZX         ZY                                  |
|                      0.2071    0.2071    0.2071                              |
+------------------------------------------------------------------------------+

[Anders Blom]

Indeed DFT is formally a zero-temperature theory in the sense that temperature doesn't really enter the equations explicitly. But as this paper attempts, you can introduce temperature via MD or similar to let it "scramble" the coordinate in the system.

The method described in the article is one option, although using AIMD (DFT) is, as the authors also note, extremely time-consuming. Given that there are decent forcefields for TiN and even Al-Ti-N (provided in QuantumATK), I would use that to reduce the calculations by maybe 1000x in time! That might also open up options to allow the lattice constant to change with temperature, which I would imagine can influence the results quite a bit, something which this paper seems to neglect.

Note that the article is extracting all quantities from a simple stress-stress analysis, not the more elaborate method used in the ElasticConstants object. This is fine for cubic materials where you only have 3 independent components of the compliance matrix, but wouldn't really work for more complex cases.

Their approach is entirely possible to run in QuantumATK with some scripting to set up the distorted structures, run the MD, and extract the data.

QuantumATK does however offer a quite different, and possible significantly better, method. Instead of averaging over all the stress data from the MD simulation, we can use the SpecialThermalDisplacement algoithm to generate a supercell structure which contains a superposition of all phonon modes at once (for different temperatures). You can then just run the ElasticConstants analysis on a single structure, and repeat this for different temperatures. The only caveat is that you need to first compute the DynamicalMatrix (i.e. the phonons) for the supercell, but with a forcefield this is very fast.

So the workflow would be

1. Generate structure with SQS (if you are looking at a random alloy, else just make a normal supercell)
2. Compute the DynamicalMatrix using a forcefield
3. Apply the SpecialThermalDisplacement
4. Compute ElasticConstants
Repeat for each temperature and/or molar composition for alloys.

I would also check the lattice constant as a function of temperature by relaxing the supercell (with the forcefield) without constraining the lattice constant.

All in all this would be a significant upgrade to the referenced paper, in that it is (much) faster and more general method (and just different, which offers a nice method comparison). So you can run their AlTiN alloy and compare the results, for a nice paper! Perhaps the only benefit of the MD technique is that it capture anharmonicity, whereas the phonons are computed in a harmonic approximation, but then again we are doing quite a few other approximations, so that probably doesn't matter too much.

[Habib]

Hi Anders,
Many thanks for the detailed response. I agree with you that it has the potential to be a good article. I will try that approach and get back to you if I encounter any serious issues.

[Habib]

Hi Anders,

I followed your method and calculated the Elastic Constants of ZrC at three different temperatures (0 K, 300 K, and 1500 K) using the provided forcefield in QuantumATK (ReaxFF_CHONSiPtZrYBaTi_2013). Since these are very fast calculations, I also computed the Dynamical Matrix for all of them. However, I got similar results for all of these three temperatures.


+------------------------------------------------------------------------------+
| Elastic Constants in GPa                                                     |
+------------------------------------------------------------------------------+
|   402.38     101.43     101.43       0.00       0.00       0.00              |
|              402.38     101.43       0.00       0.00       0.00              |
|                         402.38       0.00       0.00       0.00              |
|                                    192.39       0.00       0.00              |
|                                               192.39       0.00              |
|                                                          192.39              |
+------------------------------------------------------------------------------+

The scripts of these three cases are attached. Could you please review them and let me know where I am making mistakes?

Kind regards,
Habib