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.916334491484                   -0.69    4.2    742ms
  2   -7.921170697168       -2.32       -1.52    1.0    333ms
  3   -7.921406229071       -3.63       -2.48    1.4    343ms
  4   -7.921440911666       -4.46       -2.83    2.3    482ms
  5   -7.921441741008       -6.08       -3.12    1.0    313ms
  6   -7.921442006845       -6.58       -4.54    1.0    311ms
  7   -7.921442022075       -7.82       -5.08    2.6    534ms
  8   -7.921442022139      -10.20       -5.43    1.8    386ms
  9   -7.921442022143      -11.33       -6.05    1.1    322ms
 10   -7.921442022144      -12.36       -5.98    1.7    387ms
 11   -7.921442022144      -12.51       -7.00    1.0    315ms
 12   -7.921442022144      -13.85       -7.30    2.0    441ms
 13   -7.921442022144      -14.75       -7.97    1.0    322ms
 14   -7.921442022144      -15.05       -8.94    1.6    380ms

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.85s  Non-interacting
   1        0       1        1.03     15.2   12.8s  
   2        0       2        0.53     14.5   1.26s  
   3        0       3       -0.01     14.1   1.14s  
   4        0       4       -0.87     13.8   1.13s  
   5        0       5       -1.74     11.8   1.01s  
   6        0       6       -2.43     10.0   897ms  
   7        0       7       -3.33      8.8   805ms  
   8        0       8       -4.02      7.8   746ms  
   9        0       9       -5.04      6.4   642ms  
  10        0      10       -6.06      3.9   479ms  
  11        0      11       -7.08      2.3   385ms  
  12        0      12       -8.27      1.0   300ms  
  13        1       1       -7.53     17.8   1.82s  Restart
  14        1       2       -7.87      3.1   399ms  
  15        1       3       -8.71      2.1   349ms  
                                      19.2   2.12s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.39s  Non-interacting
   1        0       1        1.03     15.2   1.30s  
   2        0       2        0.53     14.5   1.17s  
   3        0       3       -0.01     14.1   1.14s  
   4        0       4       -0.87     13.7   1.12s  
   5        0       5       -1.74     11.8   1.01s  
   6        0       6       -2.43     10.0   881ms  
   7        0       7       -3.33      8.9   803ms  
   8        0       8       -4.02      7.8   740ms  
   9        0       9       -5.04      6.4   646ms  
  10        0      10       -6.06      3.9   473ms  
  11        0      11       -7.08      2.3   373ms  
  12        0      12       -8.27      1.0   300ms  
  13        1       1       -7.53     17.8   1.62s  Restart
  14        1       2       -7.89      3.1   383ms  
  15        1       3       -8.70      2.0   330ms  
                                      19.2   2.35s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      19.9   1.40s  Non-interacting
   1        0       1        1.03     15.2   1.28s  
   2        0       2        0.53     14.5   1.15s  
   3        0       3       -0.01     14.1   1.14s  
   4        0       4       -0.87     13.7   1.15s  
   5        0       5       -1.74     11.8   999ms  
   6        0       6       -2.43     10.0   881ms  
   7        0       7       -3.33      8.9   803ms  
   8        0       8       -4.02      7.8   757ms  
   9        0       9       -5.04      6.4   639ms  
  10        0      10       -6.06      3.9   490ms  
  11        0      11       -7.08      2.2   386ms  
  12        0      12       -8.27      1.0   298ms  
  13        1       1       -7.53     17.9   1.61s  Restart
  14        1       2       -7.87      3.1   380ms  
  15        1       3       -8.67      2.1   337ms  
                                      19.2   1.48s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.37s  Non-interacting
   1        0       1        1.03     15.2   1.35s  
   2        0       2        0.53     14.5   1.17s  
   3        0       3       -0.01     14.1   1.14s  
   4        0       4       -0.87     13.8   1.12s  
   5        0       5       -1.74     11.8   1.00s  
   6        0       6       -2.43     10.0   881ms  
   7        0       7       -3.33      8.9   805ms  
   8        0       8       -4.02      7.8   750ms  
   9        0       9       -5.04      6.4   634ms  
  10        0      10       -6.06      3.9   469ms  
  11        0      11       -7.08      2.3   378ms  
  12        0      12       -8.27      1.0   301ms  
  13        1       1       -7.52     17.8   1.60s  Restart
  14        1       2       -7.86      3.1   401ms  
  15        1       3       -8.70      2.1   1.25s  
                                      19.2   1.46s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.38s  Non-interacting
   1        0       1        1.03     15.2   1.28s  
   2        0       2        0.53     14.5   1.18s  
   3        0       3       -0.01     14.1   1.14s  
   4        0       4       -0.87     13.7   1.14s  
   5        0       5       -1.74     11.8   1.00s  
   6        0       6       -2.43     10.0   886ms  
   7        0       7       -3.33      8.9   804ms  
   8        0       8       -4.02      7.8   738ms  
   9        0       9       -5.04      6.4   634ms  
  10        0      10       -6.06      3.9   470ms  
  11        0      11       -7.08      2.2   375ms  
  12        0      12       -8.27      1.0   295ms  
  13        1       1       -7.53     17.7   1.60s  Restart
  14        1       2       -7.91      3.1   374ms  
  15        1       3       -8.70      2.0   334ms  
                                      19.2   1.48s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.37s  Non-interacting
   1        0       1        1.03     15.2   1.28s  
   2        0       2        0.53     14.5   1.17s  
   3        0       3       -0.01     14.1   1.15s  
   4        0       4       -0.87     13.7   1.12s  
   5        0       5       -1.74     11.8   998ms  
   6        0       6       -2.43     10.0   873ms  
   7        0       7       -3.33      8.9   802ms  
   8        0       8       -4.02      7.8   741ms  
   9        0       9       -5.04      6.4   636ms  
  10        0      10       -6.06      3.9   467ms  
  11        0      11       -7.08      2.2   381ms  
  12        0      12       -8.27      1.0   305ms  
  13        1       1       -7.52     17.8   1.63s  Restart
  14        1       2       -7.89      3.1   1.23s  
  15        1       3       -8.70      2.1   348ms  
                                      19.2   1.60s  Final orbitals
134.705877 seconds (101.25 M allocations: 30.776 GiB, 4.46% gc time, 25.73% 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