Noncollinear (spin_polarization=:full) calculation added#1273
Noncollinear (spin_polarization=:full) calculation added#1273physarief78 wants to merge 3 commits into
spin_polarization=:full) calculation added#1273Conversation
|
Thanks! It looks quite good! Did you test it with a genuinely non-collinear system to see if it works well? Ideally we should check the result against an existing code (abinit or quantum espresso). The main problem here is that, in the DFTK philosophy, we want to hardcode as less about DFT as possible. The motivation for this is that there are many other possible uses of the multicomponent DFT-like models : elasticity, maxwell's equations, and some coarse-grained models of graphene. A previous attempt was made in #851 but stalled. That was some time ago but we discussed a lot about data structures at the time: the constraints are that a) we don't want to penalize too much code readability (ie the code should not have a lot of The most important choices we have to make are how to represent orbitals and densities. We settled in #851 on psi[sigma, G, n] I think (but that was some time ago, before we had collinear spin). In #849 (comment) I said psi[G,n][sigma] with svectors which is memory-layout compatible but not exactly the same as psi[sigma, G, n]. Looking at the code now it looks like you've basically hardcoded the representation of psi in the functions that use it (ie you made no modification to fft.jl), but it looks like in memory layout it might be closer to psi[G, sigma, n]? We should really make it a part of the "wavefunction API". Regarding densities, it's not completely clear to me what the correct representation should be in the general case: probably it should be density matrices (a 2x2 hermitian matrix). Is it very bad to just use and compute this in the uncompressed zeros(ComplexF64, 2, 2) form? If so we should probably have helper functions that convert to and from the 2x2 hermitian format to the R^3 format. @mfherbst what do you think? |
|
Thank you for the thoughtful feedback and suggestions. To be honest, I haven't yet had the chance to test the implementation against existing codes like Quantum ESPRESSO or ABINIT. This is mainly because I've been focusing on getting Regarding the code structure, I did try to keep it as simple, clear, and general as possible. However, I encountered quite a few conflicts along the way, and at some point it became mentally taxing to balance generality with functionality. That's why, in the current PR, I ended up hardcoding certain parts—partly as a way to move forward and gather early feedback. I realize this goes against the design philosophy of DFTK. Your insights are very helpful, and I will now work on refactoring the code to better align with the generic, multicomponent framework you described. I'll also try to incorporate the ideas about wavefunction and density representations we discussed, and I'll keep the updated on progress and any remaining challenges. One aspect I'm still struggling with is the mapping between the density matrix and the magnetization vector. It's not a straightforward one-to-one mapping if we want to preserve the full Hermitian structure for physical robustness. |
|
Why not 1 to 1? Is it not just splitting a 2x2 hermitian matrix onto the four Pauli matrices (with sigma0=Id)? |
|
Ah, you are completely right, and I phrased that poorly! My apologies. Mathematically, it is absolutely a 1-to-1 mapping via the Pauli matrices When I said I was struggling with the mapping and 'robustness,' I didn't mean the mathematical theory. I meant I was struggling with the best way to implement this cleanly in Julia while adhering to DFTK's generic architecture. If we just use an uncompressed However, your suggestion makes perfect sense. About to implement the helper functions to explicitly convert between the 2 x 2 Hermitian matrix form and the |
|
That decomposition is fully general, it's just selecting a real basis for the space of traceless hermitian matrices. In the 2x2 case it's R^3 (two because complex times one for the off-diagonal, two for the diagonal, minus one for the trace), in the 3x3 case it's 8 (two times three for the off-diagonal, three for the diagonal, minus one for the trace), etc |
|
Ah, I see exactly what you mean now! I was getting tunnel vision on the Pauli matrices as a special 2 x 2 case and completely missed the broader algebraic structure. Thank you for spelling that out. So since the decomposition into a scalar trace and an |
2 participants
Dear maintainers and community
According to DFTK issue #1272, I have implemented support for noncollinear calculations (
spin_polarization=:full) in DFTK. To achieve this, as can be seen in the changed files, 21 files were modified. This was necessary because noncollinear calculations involve the x, y, and z components of the spin direction; in other words, the spin direction is no longer restricted to just up or down, but can point in any direction. Furthermore, we also need to compute the probability density, which together with the three spin components results in a four-dimensional problem. Since DFTK previously assumed a simplified definition of spin polarization (i.e., collinear up/down) before the support forspin_polarization=:fullwas added, many changes were required to enable noncollinear calculations properly and safely. "Safely" here means that these changes should not affect the existing modes such asspin_polarization=:collinear,:spinless, or:none, nor any other solvers in general.Getting more detailed we can see this mathematically by comparing the collinear and noncollinear calculation. The fundamental difference lies in how we treat the quantum mechanical spin state of the electrons.
In collinear DFT, the global quantization axis (usually the$z$ -axis) is fixed. An electron is strictly either "spin-up" or "spin-down". The wavefunctions are simple scalar fields for each spin channel:
In non-collinear DFT, the spin quantization axis can rotate freely from point to point in the crystal. Because spin is a two-state quantum system, we must promote the wavefunction to a 2-component Pauli spinor. Every electronic state$i$ becomes a mixture of up and down:
Because the wavefunctions are vectors, the electron density is no longer a simple scalar$n(\mathbf{r})$ . It becomes a $2 \times 2$ Hermitian density matrix at every point in space:
Any$2 \times 2$ Hermitian matrix can be decomposed into a scalar charge density $n(\mathbf{r})$ and a 3D magnetization vector $\mathbf{m}(\mathbf{r})$ using the Pauli matrices ($\boldsymbol{\sigma} = (\sigma_x, \sigma_y, \sigma_z)$):
This is exactly why we now store the density as a four-component array:$\rho[:,:,:,1]$ is $n$ , and $\rho[:,:,:,2:4]$ are $m_x, m_y, m_z$ .
The effective Kohn-Sham Hamiltonian dictates how the orbitals evolve:
In collinear DFT, the exchange-correlation potential$\hat{V}_{\text{xc}}$ only has a $z$ -component. The Hamiltonian is block-diagonal, meaning we essentially solve two smaller, separate equations (one for up, one for down):
In non-collinear DFT, the exchange-correlation potential acts like a local magnetic field$\mathbf{B}_{\text{xc}}(\mathbf{r})$ that points in arbitrary 3D directions. The Hamiltonian becomes a fully coupled $2 \times 2$ matrix:
Because of those off-diagonal terms ($B_x \pm iB_y$ ), an electron's spin can "flip" or rotate as it moves through the lattice. This fully coupled matrix is twice as large, which is why non-collinear calculations take significantly more RAM and CPU time to diagonalize.
Standard exchange-correlation functionals (like LDA or PBE) are only mathematically defined for scalar up/down densities, not vector fields. To use libxc, we must mathematically project the non-collinear state into a local collinear state at every single grid point. This is exactly what we implemented in our xc.jl file:
We compute the magnitude of the local magnetization vector:
Then we construct the local up/down densities by diagonalizing the density matrix:
Step B: Call libxc$n_+$ and $n_-$ into libxc to get the scalar potentials $V_+$ and $V_-$ .
We feed
Step C: Rotate the potential back to the global frame$2 \times 2$ matrix form to build the exchange-correlation magnetic field:
Using the chain rule, we map the scalar potentials back into the
This maps directly to the$V_n$ and $\text{potential}[i, 2:4] \mathrel{+}= V_m \cdot \mathbf{m} / |\mathbf{m}|$ lines in our Julia implementation. By following these steps, we ensure that the noncollinear calculation correctly handles the spin degrees of freedom while still leveraging the standard libxc library.
To check the results we can do this step-by-step:
printed scf calculation results:
The results:
If we print the code above, this will gives us the summary result of:
Conclusion:
Based on the results, the implementation of noncollinear support (
spin_polarization=:full) was successfully carried out in DFTK by modifying 21 files to handle four-component spinor wave functions and fully bound Hamiltonians, while retaining all existing functions. This implementation is mathematically rigorous, based on the Pauli spinor formalism and a local collinear approximation for the exchange-correlation function. Numerical tests on solid iron show excellent agreement: collinear and noncollinear calculations converge to the same total energy (difference < 10⁻⁸ Hartree) and identical magnetic moment magnitude (4.366 μB), with noncollinear results also providing complete vector components (Mx, My, Mz). Further visualization of spin textures, band structures, and state densities confirms physical consistency. These features are numerically stable and physically consistent.That is all for my PR. I hope this contribution will help make DFTK more robust and enable many future methodological developments. Thank you in advance for any feedback. This is the final step—I hope the PR passes all checks before being merged.
AI statements that were used for making this PR:
GitHub Copilot:
Since I didn't know where to start, I asked GitHub Copilot which files I should modify. However, after I began modifying them and implementing the calculation for
spin_polarization=:full, errors kept appearing until I had to modify all 21 files (as can be seen in the21 files changedpart on this GitHub PR feature).Gemini AI:
I used Gemini to help me code and modify the 21 files that were changed for the
spin_polarization=:fullimplementation.ChatGPT:
I used ChatGPT to summarize my derivations into compact, short conclusions to help me improve my prompting, bridging the physics theory and guiding the statements from GitHub Copilot that implemented in Gemini AI to make the necessary file changes.