stress-genericlinearalgebra.jl

# Very basic setup, useful for testing
using DoubleFloats
using GenericLinearAlgebra
using DFTK

a = 10.26
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]

model = model_atomic(lattice, atoms, symmetries=false)
kgrid = [1, 1, 1]  # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 15          # kinetic energy cutoff in Hartree
basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)

scfres = self_consistent_field(basis, tol=1e-8)

function compute_energy(scfres_ref, a)
    lattice = a / 2 * [[0 1 1.];
                       [1 0 1.];
                       [1 1 0.]]
    Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
    atoms = [Si => [ones(3)/8, -ones(3)/8]]

    model = model_atomic(lattice, atoms, symmetries=false)
    kgrid = [1, 1, 1]  # k-point grid (Regular Monkhorst-Pack grid)
    Ecut = 15           # kinetic energy cutoff in Hartree
    basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid)

    energies, H = energy_hamiltonian(basis, scfres_ref.ψ, scfres_ref.occupation; ρ=scfres_ref.ρ)
    energies.total
end

compute_energy(scfres, 10.26)

###
### Forward mode
###

using ForwardDiff
ForwardDiff.derivative(a -> compute_energy(scfres, a), 10.26)
# [ Info: Changing fft size to 32 (smallest working size for generic FFTs)
# [ Info: Changing fft size to 32 (smallest working size for generic FFTs)
# [ Info: Changing fft size to 32 (smallest working size for generic FFTs)
# ERROR: LoadError: DimensionMismatch("arrays could not be broadcast to a common size; got a dimension with lengths 27 and 32")
# Stacktrace:
#   [1] _bcs1
#     @ ./broadcast.jl:501 [inlined]
#   [2] _bcs(shape::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}}, newshape::Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}, Base.OneTo{Int64}})
#     @ Base.Broadcast ./broadcast.jl:495
#   [3] broadcast_shape
#     @ ./broadcast.jl:489 [inlined]
#   [4] combine_axes
#     @ ./broadcast.jl:484 [inlined]
#   [5] instantiate
#     @ ./broadcast.jl:266 [inlined]
#   [6] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{3}, Nothing, typeof(*), Tuple{Array{Float64, 3}, Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1}, 3}}})
#     @ Base.Broadcast ./broadcast.jl:883
#   [7] macro expansion
#     @ ~/.julia/dev/DFTK.jl/src/terms/local.jl:17 [inlined]
#   [8] ene_ops(term::DFTK.TermAtomicLocal, ψ::Vector{Matrix{ComplexF64}}, occ::Vector{Vector{Float64}}; kwargs::Base.Iterators.Pairs{Symbol, Array{Float64, 4}, Tuple{Symbol}, NamedTuple{(:ρ,), Tuple{Array{Float64, 4}}}})
#     @ DFTK ~/.julia/packages/TimerOutputs/PZq45/src/TimerOutput.jl:210
#   [9] (::DFTK.var"#111#118"{Base.Iterators.Pairs{Symbol, Array{Float64, 4}, Tuple{Symbol}, NamedTuple{(:ρ,), Tuple{Array{Float64, 4}}}}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}})(term::DFTK.TermAtomicLocal)
#     @ DFTK ./none:0
#  [10] iterate
#     @ ./generator.jl:47 [inlined]
#  [11] collect_to!(dest::Vector{NamedTuple{(:E, :ops), Tuple{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1}, Vector{DFTK.FourierMultiplication}}}}, itr::Base.Generator{Vector{Any}, DFTK.var"#111#118"{Base.Iterators.Pairs{Symbol, Array{Float64, 4}, Tuple{Symbol}, NamedTuple{(:ρ,), Tuple{Array{Float64, 4}}}}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}}}, offs::Int64, st::Int64)
#     @ Base ./array.jl:724
#  [12] collect_to_with_first!(dest::Vector{NamedTuple{(:E, :ops), Tuple{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1}, Vector{DFTK.FourierMultiplication}}}}, v1::NamedTuple{(:E, :ops), Tuple{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1}, Vector{DFTK.FourierMultiplication}}}, itr::Base.Generator{Vector{Any}, DFTK.var"#111#118"{Base.Iterators.Pairs{Symbol, Array{Float64, 4}, Tuple{Symbol}, NamedTuple{(:ρ,), Tuple{Array{Float64, 4}}}}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}}}, st::Int64)
#     @ Base ./array.jl:702
#  [13] collect(itr::Base.Generator{Vector{Any}, DFTK.var"#111#118"{Base.Iterators.Pairs{Symbol, Array{Float64, 4}, Tuple{Symbol}, NamedTuple{(:ρ,), Tuple{Array{Float64, 4}}}}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}}})
#     @ Base ./array.jl:683
#  [14] macro expansion
#     @ ~/.julia/packages/TimerOutputs/PZq45/src/TimerOutput.jl:210 [inlined]
#  [15] macro expansion
#     @ ~/.julia/dev/DFTK.jl/src/terms/Hamiltonian.jl:124 [inlined]
#  [16] energy_hamiltonian(basis::PlaneWaveBasis{ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1}}, ψ::Vector{Matrix{ComplexF64}}, occ::Vector{Vector{Float64}}; kwargs::Base.Iterators.Pairs{Symbol, Array{Float64, 4}, Tuple{Symbol}, NamedTuple{(:ρ,), Tuple{Array{Float64, 4}}}})
#     @ DFTK ~/.julia/packages/TimerOutputs/PZq45/src/TimerOutput.jl:210
#  [17] compute_energy(scfres_ref::NamedTuple{(:ham, :basis, :energies, :converged, :ρ, :eigenvalues, :occupation, :εF, :n_iter, :n_ep_extra, :ψ, :diagonalization, :stage), Tuple{Hamiltonian, PlaneWaveBasis{Float64}, Energies{Float64}, Bool, Array{Float64, 4}, Vector{Vector{Float64}}, Vector{Vector{Float64}}, Float64, Int64, Int64, Vector{Matrix{ComplexF64}}, NamedTuple{(:λ, :X, :residual_norms, :iterations, :converged, :n_matvec), Tuple{Vector{Vector{Float64}}, Vector{Matrix{ComplexF64}}, Vector{Vector{Float64}}, Vector{Int64}, Bool, Int64}}, Symbol}}, a::ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1})
#     @ Main ~/Documents/cloud/sandbox/julia/dftk/DifferentiableDFTK/sandbox/silicon-stress/stress-total.jl:32
#  [18] (::var"#9#10")(a::ForwardDiff.Dual{ForwardDiff.Tag{var"#9#10", Float64}, Float64, 1})
#     @ Main ~/Documents/cloud/sandbox/julia/dftk/DifferentiableDFTK/sandbox/silicon-stress/stress-total.jl:43
#  [19] derivative(f::var"#9#10", x::Float64)
#     @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/derivative.jl:14
#  [20] top-level scope
#     @ ~/Documents/cloud/sandbox/julia/dftk/DifferentiableDFTK/sandbox/silicon-stress/stress-total.jl:43