Phonon computations

This is a quick sketch how to run a simple phonon calculation using DFTK.

Preliminary implementation

Practical phonon computations have only seen rudimentary testing as of now. As of now we do not yet recommend relying on this feature for production calculations. Some of the limitations are:

  • symmetries must be disabled (pass symmetries=false to the model)
  • only LDA functionals are supported
  • non-linear core corrections from the pseudopotentials are not supported
  • MPI parallelization over k-points is not supported (due to $k$ and $k+q$ interactions)

We appreciate any issues, bug reports or PRs.

First we run an SCF calculation.

using AtomsBuilder
using DFTK
using Printf
using PseudoPotentialData

pseudopotentials = PseudoFamily("cp2k.nc.sr.lda.v0_1.semicore.gth")
# Make sure to disable symmetries:
model  = model_DFT(bulk(:Si); pseudopotentials, functionals=LDA(), symmetries=false)
basis  = PlaneWaveBasis(model; Ecut=10, kgrid=[4, 4, 4])
scfres = self_consistent_field(basis, tol=1e-8);
nothing  # hide
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -7.916333386731                   -0.69    4.2    712ms
  2   -7.921169437687       -2.32       -1.52    1.0    323ms
  3   -7.921407826597       -3.62       -2.48    1.5    402ms
  4   -7.921440900397       -4.48       -2.84    2.3    486ms
  5   -7.921441761122       -6.07       -3.15    1.1    322ms
  6   -7.921442007062       -6.61       -4.45    1.0    298ms
  7   -7.921442022022       -7.83       -4.91    2.5    558ms
  8   -7.921442022131       -9.96       -5.19    1.7    378ms
  9   -7.921442022142      -10.94       -5.77    1.1    339ms
 10   -7.921442022144      -11.73       -6.43    1.3    360ms
 11   -7.921442022144      -12.76       -7.32    1.9    438ms
 12   -7.921442022144      -14.15       -7.45    2.0    493ms
 13   -7.921442022144   +  -14.75       -7.99    1.0    328ms
 14   -7.921442022144   +  -14.75       -8.73    1.3    355ms

Next we compute the phonon modes at the q-point [1/4, 1/4, 1/4].

phret_q0 = @time DFTK.phonon_modes(scfres; q=[0.25, 0.25, 0.25]);
nothing  # hide
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   2.09s  Non-interacting
   1        0       1        1.03     15.2   13.6s  
   2        0       2        0.53     14.5   1.44s  
   3        0       3       -0.01     14.1   1.37s  
   4        0       4       -0.87     13.7   2.28s  
   5        0       5       -1.74     11.8   1.14s  
   6        0       6       -2.43     10.0   1.06s  
   7        0       7       -3.33      8.9   925ms  
   8        0       8       -4.02      7.8   864ms  
   9        0       9       -5.04      6.5   759ms  
  10        0      10       -6.06      3.9   569ms  
  11        0      11       -7.08      2.3   438ms  
  12        0      12       -8.27      1.0   345ms  
  13        1       1       -7.53     17.8   2.10s  Restart
  14        1       2       -7.89      3.2   461ms  
  15        1       3       -8.74      2.2   402ms  
                                      19.2   2.35s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.62s  Non-interacting
   1        0       1        1.03     15.2   1.55s  
   2        0       2        0.53     14.6   1.42s  
   3        0       3       -0.01     14.1   1.38s  
   4        0       4       -0.87     13.8   1.35s  
   5        0       5       -1.74     11.8   1.22s  
   6        0       6       -2.43     10.0   1.06s  
   7        0       7       -3.33      8.9   967ms  
   8        0       8       -4.02      7.8   877ms  
   9        0       9       -5.04      6.5   761ms  
  10        0      10       -6.06      3.9   564ms  
  11        0      11       -7.08      2.3   438ms  
  12        0      12       -8.28      1.0   349ms  
  13        1       1       -7.52     17.8   1.93s  Restart
  14        1       2       -7.89      3.2   445ms  
  15        1       3       -8.70      2.1   394ms  
                                      19.2   1.75s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.61s  Non-interacting
   1        0       1        1.03     15.2   2.47s  
   2        0       2        0.53     14.6   1.32s  
   3        0       3       -0.01     14.1   1.34s  
   4        0       4       -0.87     13.8   1.32s  
   5        0       5       -1.74     11.8   1.18s  
   6        0       6       -2.43     10.0   1.05s  
   7        0       7       -3.33      8.9   968ms  
   8        0       8       -4.02      7.9   878ms  
   9        0       9       -5.04      6.5   762ms  
  10        0      10       -6.06      3.9   566ms  
  11        0      11       -7.08      2.3   446ms  
  12        0      12       -8.27      1.0   349ms  
  13        1       1       -7.52     17.8   1.95s  Restart
  14        1       2       -7.89      3.2   438ms  
  15        1       3       -8.67      2.1   386ms  
                                      19.2   1.78s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.61s  Non-interacting
   1        0       1        1.03     15.2   1.56s  
   2        0       2        0.53     14.6   1.42s  
   3        0       3       -0.01     14.1   1.37s  
   4        0       4       -0.87     13.7   1.37s  
   5        0       5       -1.74     11.8   1.20s  
   6        0       6       -2.43     10.0   1.05s  
   7        0       7       -3.33      8.9   964ms  
   8        0       8       -4.02      7.8   886ms  
   9        0       9       -5.04      6.5   765ms  
  10        0      10       -6.06      3.9   565ms  
  11        0      11       -7.08      2.3   448ms  
  12        0      12       -8.27      1.0   351ms  
  13        1       1       -7.52     17.8   1.92s  Restart
  14        1       2       -7.88      3.2   452ms  
  15        1       3       -8.70      2.2   478ms  
                                      19.2   2.56s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.47s  Non-interacting
   1        0       1        1.03     15.2   1.49s  
   2        0       2        0.53     14.6   1.36s  
   3        0       3       -0.01     14.1   1.37s  
   4        0       4       -0.87     13.8   1.36s  
   5        0       5       -1.74     11.8   1.19s  
   6        0       6       -2.43     10.0   1.05s  
   7        0       7       -3.33      8.9   961ms  
   8        0       8       -4.02      7.8   879ms  
   9        0       9       -5.04      6.4   762ms  
  10        0      10       -6.06      3.9   553ms  
  11        0      11       -7.08      2.3   440ms  
  12        0      12       -8.28      1.0   351ms  
  13        1       1       -7.52     17.9   1.91s  Restart
  14        1       2       -7.87      3.2   456ms  
  15        1       3       -8.66      2.1   380ms  
                                      19.2   1.77s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.63s  Non-interacting
   1        0       1        1.03     15.2   1.54s  
   2        0       2        0.53     14.6   1.42s  
   3        0       3       -0.01     14.1   1.39s  
   4        0       4       -0.87     13.8   1.36s  
   5        0       5       -1.74     11.8   1.22s  
   6        0       6       -2.43     10.0   1.06s  
   7        0       7       -3.33      8.9   957ms  
   8        0       8       -4.02      7.9   897ms  
   9        0       9       -5.04      6.4   767ms  
  10        0      10       -6.06      3.9   567ms  
  11        0      11       -7.08      2.3   445ms  
  12        0      12       -8.27      1.0   349ms  
  13        1       1       -7.52     17.8   1.92s  Restart
  14        1       2       -7.89      3.2   535ms  
  15        1       3       -8.64      2.1   1.16s  
                                      19.2   1.62s  Final orbitals
153.696822 seconds (143.91 M allocations: 31.733 GiB, 4.29% gc time, 23.32% compilation time)

These are the final phonon frequencies:

for (i, ω) in enumerate(phret_q0.frequencies)
    @printf("Mode %2d: %8.3f cm-1\n", i, ω .* DFTK.hartree_to_cm⁻¹)
end
Mode  1:   93.616 cm-1
Mode  2:   93.616 cm-1
Mode  3:  223.147 cm-1
Mode  4:  478.527 cm-1
Mode  5:  484.617 cm-1
Mode  6:  484.617 cm-1