What method does GULP use to build the dynamical matrix?


I am curious how GULP calculates the dynamical matrix for phonon calculations? Does it use finite differences (direct method)? I tried looking in the documentation for this, but only found the definition of the dynamical matrix and that GULP uses ALAMODE for the third order force constants.


Hi Ethan,
The details of how the dynamical matrix is calculated are given in the manual & papers about GULP. In short, the matrix is computed analytically for most potential functions. Finite differences can be used (at gamma only) for debugging and checking, but otherwise you would only use this if analytic second derivatives are not available.

1 Like

Maybe a naive question, I’m still trying to wrap my head around anharmonic IFCs, but why are 3rd and 4th order force constants not calculated analytically? Like why use ALAMODE which fits the force constants for 3rd order?


Writing down 3rd and 4th derivatives is straight forward for a 2 atom system & so it might seem simple enough. However, when you have an N-particle system and a many-body formalism this is not at all straight forward since you end up with a 3 or 4 dimensional matrix and a lot of book keeping to do! I’d recommend trying it yourself and you’ll soon be convinced that it’s not as easy as it might sound. You might also blow the memory of your computer trying to store the fourth derivative tensor.
It should be noted that GULP does compute analytic third derivatives for a range of models (such as ionic force fields) so that you can do analytic free energy minimisation within the quasiharmonic approximation. However, the third derivatives are used on the fly to avoid the storage problem and not printed out (which would be ugly). Coupling to ALAMODE provides a convenient way of getting the thermal transport properties that is much easier than coding everything within GULP, but this could be done given a large investment of someone’s time.

Makes sense. I have already crashed my cluster a few times running out of memory for 3rd order.

I do have the impression if you have automatically differentiable code it would not matter the form of your energy potential. This would give exact derivatives regardless of the functions complexity. It is something I am looking into.

One distinction I am trying to understand though is: Would true (exact) derivatives of the energy actually re-construct the total system energy or is it only through fitting to data that you can attain an accurate re-construction.

I agree that automatic differentiation might be a useful tool to make higher level derivatives manageable. We’re also looking at JAX-MD within the group for some purposes. How well it will handle a general complex manybody force field, we’ll see…

The exact derivatives should always be more faithful to the energy surface than fitted data, with the possible exception of a very noisy, discontinuous surface where fitting might smooth out the underlying useful surface.