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

[Anders Blom]

In your scripts, the elastic constants are computed using "bulk_configuration", but the thermally excited geometry is "new_configuration".

[Habib]

Hi Anders,

Thanks for the corrections, but still the results are not convincing, i am trying my best let see.

I have another idea regarding the simulation of elastic constants at different temperatures using a combination of DFT and MD simulations (maybe i am completely wrong). 
How about this approach: 
1. First, optimize the model using DFT. 
2. Then, use the DFT-optimized model for MD simulations at the desired temperature (e.g., heat up to 1500 K or any Temperature) with the NPT ensemble. 
3. Finally, take the last snapshot of the MD-simulated model (at the desired temperature) and calculate the elastic constants using DFT again, without further optimization. 

[Anders Blom]

The concern I have with that is that whichever MD image you pick, it will have a lot of randomness to it. That is, picking 2 or 3 different images (spaced out decently in time, but in the steady-state phase) you might get quite different results. So it is probably necessary to sample several images, which not only gives you an average but also a variation, i.e. error bars.

The idea with the thermal displacement method is that it provides a single structure which is a statistical ensemble for a given temperature and thus eliminates the need for the averaging (that is, you can do just a single DFT calculation for the elastic constants). The "price" you pay is a more expensive calculation of the dynamical matrix, but with a forcefield it's not an issue anyway.

Whether or not these two methods give the same result would of course be interesting to compare! Formally they are supposed to...

Finally, keep in mind that the accuracy of the forcefield should perhaps be checked for the specific material. I doubt ZrC was specifically included when they derived the ReaxFF potenial... (In fact, the original paper is for SiC, and I think all those metal elements are just included as possible impurities.) This might be an example where no forcefield exists and you have to resort to DFT. Or, why not try the MACE or M3GNet neural network potential (included in W-2024.09) - that might actually work really well here! You can also fit a machine-learned forcefield (MTP) using QuantumATK - it should actually be relatively easy for ZrC if you "only" need phonons and not melted or amorphous structures.

Also, when using the displacement method you need to use a larger unit cell (and you also need that for the MD method). For anything involving finite temperature, you cannot operate on small unit cells; a simple way to understand that is that the phonon wavelengths are typically much larger than lattice constants (at least a few nanometers, typically).

[Habib]

Hi Anders,

Many thanks for your time.
Yes, I completely agree that the approach of selecting several MD images, calculating their average, and then comparing them is very time-consuming.
This time, I extended the model size to approximately 6 nm along the c-direction and repeated the calculations at 0 K, 300 K, 1500K, and 2000K using the same method. However, I got similar results for all different Temperatures. The scripts are attached. I think this could be due to the limitations of the forcefield. I am further increasing the size of the system to give it another try; let's see.
+------------------------------------------------------------------------------+
| Elastic Constants in GPa                                                     |
+------------------------------------------------------------------------------+
|   363.46      35.09      35.07       0.01      -0.00      -0.01              |
|              363.46      35.08       0.00      -0.02      -0.02              |
|                         363.46      -0.01      -0.01       0.01              |
|                                     50.53      -0.01      -0.01              |
|                                                50.52       0.01              |
|                                                           50.52              |
+------------------------------------------------------------------------------+

+------------------------------------------------------------------------------+
| Material properties calculated from the elastic constants:                   |
+------------------------------------------------------------------------------+
| Moduli in units of GPa:                                                      |
+------------------------------------------------------------------------------+
|                 Reuss     Voigt     Hill                                     |
+------------------------------------------------------------------------------+
| Bulk modulus:   144.5405  144.5405  144.5405                                 |
| Shear modulus:   69.8703   95.9892   82.9298                                 |
+------------------------------------------------------------------------------+
|                     X         Y         Z                                    |
| Young's modulus:   357.2838  357.2880  357.2832                              |
+------------------------------------------------------------------------------+

It seems the current forcefield calculations may not be well-suited for my system, as my system is composed of high-entropy alloys (Ti, Zr, Nb, Hf, Ta, etc).
On the other hand, simulating the dynamical matrix with DFT is very time-consuming, even for small systems. Is there any way to compute a reasonable dynamical matrix using DFT with reduced computational time?
Regarding the use of MACE or M3GNet neural network potentials, I currently don’t have access to them because I am using the 2019 version of QuantumATK. However, we have ordered the latest version through the Microelectronics Centre (STFC), although the process is taking some time. Once I gain access to the updated version, I will certainly explore these methods as well.

[Anders Blom]

Is there any way to compute a reasonable dynamical matrix using DFT with reduced computational time?

Not really. Sure, you can choose a minimal basis set (SingleZeta, even) but once you make temperature-dependent configurations, you will have 100s of atoms (and yes, you can use a Wigner-Seitz cell then, or simply sample the Gamma point only). But it will still be very time-consuming.

I think your best bet will be the neural network potentials. Of course they need to be checked for accuracy for the systems in question, but I would generally say that thermal studies is precisely one of the exciting new application areas that modern machine-learned forcefield open up!