Add NLCC support for phonon computations#1311
Conversation
|
I had to remove |
antoine-levitt
left a comment
antoine-levitt
left a comment
There was a problem hiding this comment.
Haven't checked everything but looks reasonable!
| in reduced coordinates: | ||
| ``` | ||
| dynmat[β, t, α, s] = ∂²E/∂u_sα(-q)∂u_tβ(q). | ||
| ``` |
There was a problem hiding this comment.
really not a fan of this interpretation, defining it as the derivative of the R -> F map seems much simpler to me...
There was a problem hiding this comment.
That's basically what it is since dE/du is the force. But with an extra q dependency here. The main reason I put the formula there is to know in which order to write the dynmat entries, otherwise you can get the transpose by mistake.
There was a problem hiding this comment.
right but we don't use this notation elsewhere so maybe deltaF[alpha, s] = dynmat[] deltaR[]...? Any particular reason why it's beta alpha rather than alpha beta?
There was a problem hiding this comment.
I put it in that order such that the RHS matches Baroni, and it also makes sense to me with Julia's column-major indexing to write it like that, no?
There was a problem hiding this comment.
I mean it's usually xi = Aij xj, not xj = Aji xi
| for s = 1:n_atoms, α = 1:basis.model.n_dim | ||
| # Get δH ψ | ||
| δHψs_αs = compute_δHψ_αs(basis, ψ, α, s, q) | ||
| δHψs_αs = compute_δHψ_αs(basis, ψ, α, s, q; ρ) |
There was a problem hiding this comment.
We really need this state refactor business...
| # however since the potential is real, we conjugate to obtain the -q term: | ||
| # δV = Kxc (∂ρ/∂u_sα(-q) + ∂ρcore/∂u_sα(-q)) | ||
| # = conj( Kxc (∂ρ/∂u_sα( q) + ∂ρcore/∂u_sα( q)) ) | ||
| δV_αs = conj.(apply_kernel(term, basis, δρs[α, s] + δρcores[α, s]; ρ, q)) |
There was a problem hiding this comment.
Haven't followed the complete derivation but I'm scared by this line; the whole point of the explanation I gave in #1310 is that we never conjugate q stuff. Are you sure about this conj?
There was a problem hiding this comment.
Yes 100% sure. It's also there in the local term, just hidden by Parseval's theorem. For a real function (the density is real) the conjugate should always be fine. You only get problems if you directly conjugate the psis since those are complex. I am not sure how DFPT works for spin orbit coupling or magnetic fields. Surely Dal Corso has some paper about it, but I haven't read it.
|
@antoine-levitt do you want to review this or should I ask Michael? :D |
antoine-levitt
left a comment
antoine-levitt
left a comment
There was a problem hiding this comment.
LGTM modulo the comments
| in reduced coordinates: | ||
| ``` | ||
| dynmat[β, t, α, s] = ∂²E/∂u_sα(-q)∂u_tβ(q). | ||
| ``` |
There was a problem hiding this comment.
right but we don't use this notation elsewhere so maybe deltaF[alpha, s] = dynmat[] deltaR[]...? Any particular reason why it's beta alpha rather than alpha beta?
| # ∂E/∂u_tβ(q) = ∫ Vxc ∂ρcore/∂u_tβ(q) | ||
| # Differentiate again wrt. u_sα(-q) to get the dynamical matrix element: | ||
| # ∂²E/∂u_sα(-q)∂u_tβ(q) = ∫ Kxc (∂ρ/∂u_sα(-q) + ∂ρcore/∂u_sα(-q)) (∂ρcore/∂u_tβ(q)) | ||
| # + ∫ Vxc ∂²ρcore/∂u_sα(-q)∂u_tβ(q) |
There was a problem hiding this comment.
I'm still not very comfortable with this d2E business, and it contrasts with the way everything else is written (as dR -> dF). Can't you write it like that?
There was a problem hiding this comment.
This is how Baroni writes it. I don't like to just use dF and dR as it doesn't make it clear which has q and which has -q. This derivation also directly matches the implementation
There was a problem hiding this comment.
What I don't understand is what's up with the -q? The supercell mapping from dR to dF is cell-periodic, and therefore bloch theorem applies so a displacement dR_q e^iqx produces a force dF_q e^iqx, and we compute the mapping dR_q -> dF_q. No need to involve -q. Should we zoom to discuss?
| function compute_δHψ_αs(term::TermXc, basis::PlaneWaveBasis, ψ, α, s, q; ρ) | ||
| isnothing(term.ρcore) && return nothing | ||
|
|
||
| # With an NLCC, an atom displacement triggers a change in the XC potential. |
There was a problem hiding this comment.
that comment is not super helpful, just remove?
| # Pre-allocation of large arrays for GPU efficiency | ||
| Gs = G_vectors(basis) | ||
| ρ = to_device(basis.architecture, zeros(Complex{T}, length(Gs))) | ||
| Gqs = map(G -> G+q, G_vectors(basis)) |
There was a problem hiding this comment.
I think we call those Gpk or something somewhere else?
2 participants
Cleaner than I expected. Nice framework that you built @epolack and @antoine-levitt, especially the automated testing!