Add NLCC support for phonon computations by Technici4n · Pull Request #1311 · JuliaMolSim/DFTK.jl

Technici4n · 2026-04-26T11:15:55Z

Cleaner than I expected. Nice framework that you built @epolack and @antoine-levitt, especially the automated testing!

Technici4n · 2026-04-27T07:27:24Z

I had to remove @inferred from a test unfortunately, but I just ran the full GPU (CUDA) test suite and it still works, so that shouldn't be a problem.

antoine-levitt

Haven't checked everything but looks reasonable!

antoine-levitt · 2026-04-27T16:28:07Z

+in reduced coordinates:
+```
+dynmat[β, t, α, s] = ∂²E/∂u_sα(-q)∂u_tβ(q).
+```


really not a fan of this interpretation, defining it as the derivative of the R -> F map seems much simpler to me...

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.

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?

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?

I mean it's usually xi = Aij xj, not xj = Aji xi

I tried to write it like deltaF[alpha, s] = dynmat[] deltaR[] as well, but I had to give up because it's really confusing and unusual (it's difficult to keep track of what is R-periodic and what is not).

antoine-levitt · 2026-04-27T16:34:58Z

+        # 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))


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?

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.

Technici4n · 2026-05-18T07:32:51Z

@antoine-levitt do you want to review this or should I ask Michael? :D

antoine-levitt

LGTM modulo the comments

antoine-levitt · 2026-05-18T14:28:17Z

+in reduced coordinates:
+```
+dynmat[β, t, α, s] = ∂²E/∂u_sα(-q)∂u_tβ(q).
+```


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?

antoine-levitt · 2026-05-18T14:42:13Z

+    # ∂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)


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?

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

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?

Sure, since a displacement dR_q e^iqR produces a force dF_q e^iqR, you apply a displacement dR_q e^iqR, compute the corresponding force, and then multiply by e^-iqR to recover dF_q. So what you suggest is computing the force "normally" with d/du and then multiplying by e^-iqR. What I suggest is directly computing the force with the phase removed with d/du(-q). Does this make sense? The two should be equivalent, but I prefer the latter approach since I find the required code changes more obvious that way.

If you want we can have a call, maybe tomorrow somewhere between 13:30 and 16:00, or on Friday somewhere between 13:30 and 16:00?

antoine-levitt · 2026-05-18T14:48:06Z

    # 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))


I think we call those Gpk or something somewhere else?

The naming is not super consistent. I used Gqs because that is what the local term (my inspiration) uses.

mfherbst · 2026-06-03T09:14:21Z

    @testset "Atomic $(string(typeof(method)))" for (method, elements) in zip(methods, elements)
        basis = build_basis(elements, :none)
-        ρ = @inferred guess_density(basis, method)
+        ρ = guess_density(basis, method)


Just a remark that we had real trouble with type instabilities here, so this is why the @inferred tests exist. So this was not a GPU thing but a general performance thing. I think it's a bad sign you need to remove it again and we should think very carefully if this is the right thing to do (instead of fixing the issue)

It was added in this PR, so it's presumably safe to remove as long as DFTK still runs on the GPU: #850.

A proper fix is to make the q==0 case encoded in the type system (i.e. with a different type), which I would consider out of scope for this PR.

Well, probably (did you test ?).

In any case at least make sure we document this clearly, because with this sort of stuff you never know if does not end up biting us back some day. The density and guess densities go quite deep into the code and not having the type inferrable could easily have nasty performance regressions. So please add a comment about this here (i.e. that we had the @inferred here but had to remove it because the q handling breaks type inferability) and in the atomic superposition functions themselves with the proposed solution to have a ZeroVec3{T} type.

Hmm yeah making all uses of ρ non-inferrable is quite bad. Maybe we have to bite the bullet and actually encode q==0 in the type system

Technici4n · 2026-06-10T17:25:39Z

FWIW it would be good if this PR would be merged rather soon as I have follow-up PRs that depend on this one.

Technici4n added 4 commits April 24, 2026 11:20

Move xc phonon tests to separate file

981b0c2

Add NLCC support for phonon computations

a3b8560

Up

b639998

Remove @inferred in the test

f61384e

Remove more @inferred

98d788f

antoine-levitt reviewed Apr 27, 2026

View reviewed changes

Remove remark about commensurate q

df686fc

antoine-levitt reviewed May 18, 2026

View reviewed changes

mfherbst reviewed Jun 3, 2026

View reviewed changes

Remove comment

1330868

Uh oh!

Conversation

Technici4n commented Apr 26, 2026

Uh oh!

Technici4n commented Apr 27, 2026

Uh oh!

antoine-levitt left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Technici4n commented May 18, 2026

Uh oh!

antoine-levitt left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Technici4n commented Jun 10, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants