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.916348568606                   -0.69    4.2    695ms
  2   -7.921176700582       -2.32       -1.52    1.0    287ms
  3   -7.921409863332       -3.63       -2.48    1.5    340ms
  4   -7.921440951676       -4.51       -2.85    2.3    441ms
  5   -7.921441786697       -6.08       -3.19    1.0    301ms
  6   -7.921442008555       -6.65       -4.54    1.0    302ms
  7   -7.921442022075       -7.87       -5.05    2.7    592ms
  8   -7.921442022138      -10.20       -5.35    1.8    1.16s
  9   -7.921442022143      -11.31       -5.90    1.0    297ms
 10   -7.921442022144      -11.97       -6.13    1.3    318ms
 11   -7.921442022144      -12.66       -7.23    1.0    275ms
 12   -7.921442022144      -13.80       -7.27    2.2    476ms
 13   -7.921442022144      -14.75       -7.92    1.0    325ms
 14   -7.921442022144   +    -Inf       -9.01    1.3    307ms

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   1.88s  Non-interacting
   1        0       1        1.03     15.2   13.2s  
   2        0       2        0.53     14.6   1.32s  
   3        0       3       -0.01     14.1   1.26s  
   4        0       4       -0.87     13.7   1.22s  
   5        0       5       -1.74     11.8   1.12s  
   6        0       6       -2.43     10.0   957ms  
   7        0       7       -3.33      8.9   882ms  
   8        0       8       -4.02      7.8   815ms  
   9        0       9       -5.04      6.4   695ms  
  10        0      10       -6.06      3.9   522ms  
  11        0      11       -7.08      2.3   420ms  
  12        0      12       -8.27      1.0   329ms  
  13        1       1       -7.53     17.8   1.94s  Restart
  14        1       2       -7.90      3.2   404ms  
  15        1       3       -8.73      2.1   380ms  
                                      19.2   2.22s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.45s  Non-interacting
   1        0       1        1.03     15.2   1.39s  
   2        0       2        0.53     14.6   1.27s  
   3        0       3       -0.01     14.1   1.27s  
   4        0       4       -0.87     13.8   1.24s  
   5        0       5       -1.74     11.8   1.12s  
   6        0       6       -2.43     10.0   1.07s  
   7        0       7       -3.33      8.9   1.69s  
   8        0       8       -4.02      7.8   770ms  
   9        0       9       -5.04      6.4   670ms  
  10        0      10       -6.06      3.9   492ms  
  11        0      11       -7.08      2.3   381ms  
  12        0      12       -8.27      1.0   297ms  
  13        1       1       -7.52     17.8   1.73s  Restart
  14        1       2       -7.89      3.2   395ms  
  15        1       3       -8.71      2.1   352ms  
                                      19.2   1.62s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.50s  Non-interacting
   1        0       1        1.03     15.2   1.43s  
   2        0       2        0.53     14.6   1.15s  
   3        0       3       -0.01     14.1   1.11s  
   4        0       4       -0.87     13.8   1.10s  
   5        0       5       -1.74     11.8   978ms  
   6        0       6       -2.43     10.0   865ms  
   7        0       7       -3.33      8.9   780ms  
   8        0       8       -4.02      7.8   719ms  
   9        0       9       -5.04      6.4   618ms  
  10        0      10       -6.06      3.9   461ms  
  11        0      11       -7.08      2.3   364ms  
  12        0      12       -8.26      1.0   290ms  
  13        1       1       -7.52     17.9   1.54s  Restart
  14        1       2       -7.88      3.2   363ms  
  15        1       3       -8.67      2.1   328ms  
                                      19.2   1.44s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.33s  Non-interacting
   1        0       1        1.03     15.2   1.23s  
   2        0       2        0.53     14.6   1.13s  
   3        0       3       -0.01     14.1   1.10s  
   4        0       4       -0.87     13.7   1.08s  
   5        0       5       -1.74     11.8   961ms  
   6        0       6       -2.43     10.0   853ms  
   7        0       7       -3.33      8.9   772ms  
   8        0       8       -4.02      7.8   713ms  
   9        0       9       -5.04      6.4   609ms  
  10        0      10       -6.06      3.9   449ms  
  11        0      11       -7.08      2.3   451ms  
  12        0      12       -8.27      1.0   1.07s  
  13        1       1       -7.53     17.8   1.69s  Restart
  14        1       2       -7.89      3.2   361ms  
  15        1       3       -8.74      2.1   351ms  
                                      19.2   1.61s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.50s  Non-interacting
   1        0       1        1.03     15.2   1.42s  
   2        0       2        0.53     14.6   1.29s  
   3        0       3       -0.01     14.1   1.25s  
   4        0       4       -0.87     13.8   1.23s  
   5        0       5       -1.74     11.8   1.09s  
   6        0       6       -2.43     10.0   981ms  
   7        0       7       -3.33      8.9   890ms  
   8        0       8       -4.02      7.8   821ms  
   9        0       9       -5.04      6.4   711ms  
  10        0      10       -6.06      3.9   540ms  
  11        0      11       -7.08      2.3   418ms  
  12        0      12       -8.27      1.0   332ms  
  13        1       1       -7.53     17.7   1.77s  Restart
  14        1       2       -7.90      3.2   424ms  
  15        1       3       -8.72      2.1   358ms  
                                      19.2   1.63s  Final orbitals
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      20.0   1.48s  Non-interacting
   1        0       1        1.03     15.2   1.41s  
   2        0       2        0.53     14.6   1.29s  
   3        0       3       -0.01     14.1   1.26s  
   4        0       4       -0.87     13.8   1.23s  
   5        0       5       -1.74     11.8   1.11s  
   6        0       6       -2.43     10.0   973ms  
   7        0       7       -3.33      8.9   894ms  
   8        0       8       -4.02      7.8   820ms  
   9        0       9       -5.04      6.4   696ms  
  10        0      10       -6.06      3.9   529ms  
  11        0      11       -7.08      2.3   414ms  
  12        0      12       -8.26      1.0   333ms  
  13        1       1       -7.53     17.8   2.59s  Restart
  14        1       2       -7.88      3.2   424ms  
  15        1       3       -8.67      2.1   315ms  
                                      19.2   1.48s  Final orbitals
139.000829 seconds (101.99 M allocations: 30.814 GiB, 3.72% gc time, 25.42% 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