Automate elastic_constants by AD-DFPT as a postprocessing

[niklasschmitz]

Following up on the ideas at #1162.

This PR aims to automate the clamped-ion elastic constants as a convenient postprocessing function elastic_constants(scfres) for 1) the generic case of not using symmetry (most expensive, 6xDFPT on a fully unfolded kgrid...), 2.1) the most common case of cubic crystals (least expensive, 1xDFPT on a partially unfolded kgrid). Optimizations for other the point groups could be implemented later following a similar pattern.

progress:

• Implement generic case
• Implement cubic case
• use SCF results from scfres and partially unfold (Warm-start strained SCF from unfolded scfres in elastic_tensor #1308; Update: just warm-starting with the correct density and tightened diagonalization tolerance)
• more testing
• mention in elastic constants docs example

[mfherbst]

Something that bothers me as well is that it's basically impossible to send any kind of parameters to the DFPT solver in your implementation. There is the ResponseOptions for this purpose. Maybe just taking one as a keyword argument from the elastic constant wrapper function would be simple solution ?

[mfherbst]

Also in

strained_lattice = DFTK.voigt_strain_to_full(0.01 * strain_pattern) * model0.lattice

The displacement 0.01 should depend on the SYMMETRY_TOLERANCE

[niklasschmitz]

Something that bothers me as well is that it's basically impossible to send any kind of parameters to the DFPT solver in your implementation. There is the ResponseOptions for this purpose. Maybe just taking one as a keyword argument from the elastic constant wrapper function would be simple solution ?

Agreed. I now added the more general solution to pass arbitrary kwargs_scf... to elastic_tensor.

The displacement 0.01 should depend on the SYMMETRY_TOLERANCE

Makese sense. I added tol_symmetry as keyword argument and displace by 100 * tol_symmetry (the 1e-5 default SYMMETRY_TOLERANCE gets 1e-3 strain).

[mfherbst]

This is code duplication and I'm not sure if the strings you get out of solver and eigensolver are very useful. E.g. eigensolver could just be a large bunch of numbers (which could blow up the size of the json). Have you checked that ?

[niklasschmitz]

It seems to work fine for the function arguments, it just prints the name. On examples/silicon.jl it looks like this:

julia> dict = Dict()
Dict{Any, Any}()

julia> for sym in (:is_converged, :nbandsalg, :fermialg, :diagtolalg, :solver, :eigensolver)
     if hasproperty(scfres, sym)
         dict[string(sym)] = string(getproperty(scfres, sym))
     end
 end

julia> dict
Dict{Any, Any} with 6 entries:
  "nbandsalg"    => "AdaptiveBands(4, 7, 1.0e-6, 0.01)"
  "eigensolver"  => "lobpcg_hyper"
  "is_converged" => "ScfConvergenceDensity(1.0e-8)"
  "diagtolalg"   => "AdaptiveDiagtol(0.2, nothing, 0.005, 0.03)"
  "solver"       => "anderson"
  "fermialg"     => "FermiTwoStage()"