Hi Sherif, I just wanted to follow up to clarify a few things @mkhorton said. For future reference, please keep your questions concerning a single topic in one thread, or provide a reference to prior threads with related content. Since we have multiple members of our team answering questions, we got confused here trying to respond because it wasn’t clear from the context how you’d gotten to this point.
Although pymatgen doesn’t directly calculate debye temperature from bulk and shear modulus the way it’s implemented right now, it is doing so using only the VRH moduli derived from the elastic tensor, so you can crudely estimate the debye temperature using only these two quantities (as we’ve discussed before).
I can’t verify whether your mathematica code is quantitatively correct, but it looks like the basic structure (in which you find the speed of sound first and estimate debye temperature from that) is right. If you want to check it, you might try with more well-known material (say Fe) to see if your method reasonably approximates the experimental value or the one you can get with pymatgen and the API.
Lastly, I want to emphasize what @mkhorton said: this is a relatively crude estimate because the values of K_VRH and G_VRH are estimated based on the machine learning algorithm he linked. Also, that machine learning algorithm is known to do somewhat poorly on F-block elements like Gd, which is why a warning appears on both the website and the return data.