Phonon +q instead of -q?

[Technici4n]

The dynamical matrix is supposed to be a mixed derivative wrt one monochromatic displacement at $+q$ and one monochromatic displacement at $-q$. From the code it seems to me that we compute the $\delta \rho$ at momentum $+q$. So for example I would expect something like $\delta V_{-q} \delta \rho_{+q}$ for the local term, but instead I see that the potential change is also computed at +q?

DFTK.jl/src/terms/local.jl

Lines 193 to 197 in b6b1f7d

	# dynmat_δH, which is ∫δρδV.
	dynmat_δH = zeros(S, 3, n_atoms, 3, n_atoms)
	for s = 1:n_atoms, α = 1:n_dim
	dynmat_δH[:, :, α, s] .-= stack(forces_local(S, basis, δρs[α, s], q))
	end

Maybe we are lucky by chance and we are actually computing $D(-q)$ consistently by accident (I think this could be the case but I am not 100% sure).

[antoine-levitt]

Yeah I remember this was very confusing. I think we didn't want to take into account time reversal symmetry explicitly so we wanted to define the dynamical matrix as mapping dR to dF under a complex displacement dR. But then that mixes q and -q anyway because of the psibar. In the end there's a trick that they use in Baroni and I don't remember if we use the same one. Maybe @epolack remembers? We definitely should have written this down at the time...

[antoine-levitt]

Right the trick was that the delta rho is psibar dpsi + dpsibar psi. We don't like the second one because it couples the +q and -q phonon. So what we do is to compute it as 2 psibar dpsi. I think this is in fact fine because of a funny complex analysis reason. Surely we must have written this down somewhere...

[antoine-levitt]

So the ground truth is the supercell. There, the map χ0 = δV -> δψ -> δρ is a priori a translation invariant map from real to real. The map δV -> δψ is complex-linear, but not δψ->δρ (because δρ = ψbar δψ + δψbar ψ). So strictly speaking the Bloch theorem does not apply (because it requires complex vector spaces). We can't a priori just take δV = δVq e^iqx to compute δρ = δρq e^iqx: if we do that, we will have to compute δψ at +q and -q. However (plus or minus signs, factors and multiple bands),

χ0 = ∑_k,q' |ψkbar ψk+q'><ψkbar ψk+q'| / (εk-εk+q').

So this actually extends χ0 to a C-linear map (there's nothing preventing me to take a complex δV in this formula). If we hit it with δV = δVq e^iqx, this excites only the q'=q term, and

χ0 δVq e^iqx = ∑_k |ψkbar ψk+q><ψkbar ψk+q| / (εk-εk+q)
= ∑_k ψkbar δψk

where

δψk = ψk+q <ψkbar ψk+q| δV> / (εk-εk+q)
= e^ik+q uk+q <uk+q| δVq uk> / (εk-εk+q)
= e^ik+q δuk

where δuk solves

(Hk+q -εk) δuk = - δVq uk (up to projections, occupations, blabla)

[I think Baroni might define δψk as δψk+q, or maybe we do in the code, I never remember: the point is that δψk is a bloch wave at momentum k+q, so it's confusing whatever choice you make]

So at least that's why in the code we use δρ = ψkbar δψk, without conjugates. Does this resolve your +q / -q objection or is it something else?

[Technici4n]

I think the computation of $\delta \rho$ is correct, i.e. we sort of double the +q part to avoid having to deal with the -q components. I wonder how that works when the time reversal symmetry is broken, but that's not really a concern now I guess. And overall very confusing, totally agreed! (Especially the fact that the reponse of the k momentum is at k+q...)

My concern is about how $\delta \rho$ gets used down the line to assemble the dynamical matrix. I expect a mixed derivative wrt. +q and -q. What we can also do is take the complex conjugate of $\delta \rho$ at q to get the -q component (always true since the density is real). This is what Baroni does in Eq. (86). What confuses me is that I don't see any -q in the dynmat, nor a conjugation of the density response. Let me know if this is not clear ^^

[antoine-levitt]

I wonder how that works when the time reversal symmetry is broken

The derivation I gave above does not depend on TRS. We don't follow exactly Baroni in part because of this.

My concern is about how δρ gets used down the line to assemble the dynamical matrix. I expect a mixed derivative wrt. +q and -q.

Why do you expect this? The DM is computed from δR -> δV -> δρ -> δF (with also a direct path δR -> δF for eg Ewald). The only place that potentially mixes +q and -q is δV->δρ, and we use the trick above to not actually mix them.

[Technici4n]

The definition of the dynamical matrix should be

$$D_{s\alpha,t\beta}(q) = \frac{\partial^2 E}{\partial u_{s\alpha}^*(q)\partial u_{t\beta}(q)} = \frac{\partial^2 E}{\partial u_{s\alpha}(-q)\partial u_{t\beta}(q)}.$$

(If you think about it, we really need a q and a -q otherwise the dynamical matrix would carry some a 2q momentum)

So if we take a local potential only i.e. $E = \int V(r) \rho(r) \mathrm d r$, we get by first differentiating wrt. $u_{t\beta}(q)$ with Hellmann-Feynman and then taking the derivative again wrt. $u_{s\alpha}^*(q)$:

$$\int \frac{\partial^2 V}{\partial u_{s\alpha}^*(q)\partial u_{t\beta}(q)} \rho \mathrm d r + \int \frac{\partial V}{\partial u_{t\beta}(q)} \frac{\partial \rho}{\partial u_{s\alpha}^*(q)} \mathrm d r.$$

The second term is equal to this if we follow Baroni:

$$\int \frac{\partial V}{\partial u_{t\beta}(q)} \left(\frac{\partial \rho}{\partial u_{s\alpha}(q)}\right)^* \mathrm d r.$$

Alternatively we could first differentiate wrt. $u_{s\alpha}^*(q) = u_{s\alpha}(-q)$ and then take the second derivative wrt. $u_{t\beta}(q)$ to obtain for the second term:

$$\int \frac{\partial V}{\partial u_{s\alpha}(-q)} \frac{\partial \rho}{\partial u_{t\beta}(q)} \mathrm d r.$$

But it seems that DFTK is doing neither of the two?

[Technici4n]

Hmm I think I get it.

Use Parseval's identity:

$$\int \underbrace{\frac{\partial V}{\partial u_{t\beta}(q)}}_{\delta V} \left(\underbrace{\frac{\partial \rho}{\partial u_{s\alpha}(q)}}_{\delta \rho}\right)^* \mathrm d r = \int \delta V \delta \rho^* \mathrm d r = \sum_G \widehat{\delta V}(G) \left(\widehat{\delta \rho}(G)\right)^*.$$

This latter thing is indeed what DFTK is computing!

[Technici4n]

Thanks! ^^

[antoine-levitt]

I don't follow your derivation, for me the dynamical matrix is the q fiber of the dR to dF operator. I'm leaving the issue open as a documentation issue, I'll try to do something

[antoine-levitt]

Thinking about it more the formula I gave for chi0 relies on trs I think. Let me think about it some more...

[antoine-levitt]

Yeah Claude (opus, sonnet doesn't see the problem) says without trs we need two sternheimer and points to https://arxiv.org/pdf/1906.11673. We were not careful enough with tests it seems (or phonons predate magnetic fields). We need to add an assert and explain this better