Spin-orbit coupling and non-collinear calculations
Spin-direction constrained DFT
Suppose we want to constrain the direction of the spin magnetic moment \(\vec{m}_a\) at some atomic site \(a\) along some target direction defined by the unit vector \(\vec{\hat{u}}_a\). This can be done by introducing the penalty functional
\[E_\mathrm{cDFT}=\Lambda
\left({\vphantom{\frac{1}{1}}\vec{m}_a-\vec{\hat{u}}_a\left({\vec{\hat{u}}_a\cdot\vec{m}_a}\right)}\right)^2
=\Lambda\left({\vphantom{\frac{1}{1}}\vec{m}_a\cdot\vec{m}_a-\left({\vec{\hat{u}}_a\cdot\vec{m}_a}\right)^2}\right),\]
with the penalty \(\Lambda\) in units of eV pr. \(\mu_\mathrm{B}^2\). We take advantage of the PAW formalism to define the local magnetic moment at site \(a\) which we want to constrain. Reminder, the atomic spin density matrices are defined through
\[\begin{split}D^a_{x,i_1i_2}&=\sum_{\vec{k}n}f_{\vec{k}n}\left({\braket{\widetilde{\psi}_{\vec{k}\uparrow n}}{\widetilde{p}^a_{i_1}}\braket{\widetilde{p}^a_{i_2}}{\widetilde{\psi}_{\vec{k}\downarrow n}}+\braket{\widetilde{\psi}_{\vec{k}\downarrow n}}{\widetilde{p}^a_{i_1}}\braket{\widetilde{p}^a_{i_2}}{\widetilde{\psi}_{\vec{k}\uparrow n}}}\right), \\
D^a_{y,i_1i_2}&=-i\sum_{\vec{k}n}f_{\vec{k}n}\left({\braket{\widetilde{\psi}_{\vec{k}\uparrow n}}{\widetilde{p}^a_{i_1}}\braket{\widetilde{p}^a_{i_2}}{\widetilde{\psi}_{\vec{k}\downarrow n}}-\braket{\widetilde{\psi}_{\vec{k}\downarrow n}}{\widetilde{p}^a_{i_1}}\braket{\widetilde{p}^a_{i_2}}{\widetilde{\psi}_{\vec{k}\uparrow n}}}\right), \\
D^a_{z,i_1i_2}&=\sum_{\vec{k}n}f_{\vec{k}n}\left({\braket{\widetilde{\psi}_{\vec{k}\uparrow n}}{\widetilde{p}^a_{i_1}}\braket{\widetilde{p}^a_{i_2}}{\widetilde{\psi}_{\vec{k}\uparrow n}}-\braket{\widetilde{\psi}_{\vec{k}\downarrow n}}{\widetilde{p}^a_{i_1}}\braket{\widetilde{p}^a_{i_2}}{\widetilde{\psi}_{\vec{k}\downarrow n}}}\right).\end{split}\]
If we only wish to constrain the PAW part of the magnetic moment at site \(a\), we can write (neglecting the negative sign as is standard in GPAW)
\[m_{a,x}=\mu_\mathrm{B}\sum_{i_1i_2}D^a_{x,i_1i_2}N_{i_1i_2},\quad
m_{a,y}=\mu_\mathrm{B}\sum_{i_1i_2}D^a_{y,i_1i_2}N_{i_1i_2},\quad
m_{a,z}=\mu_\mathrm{B}\sum_{i_1i_2}D^a_{z,i_1i_2}N_{i_1i_2},\]
with
\[N_{i_1i_2}=
\braket{\phi^a_{i_1}}{\phi^a_{i_2}}_{\hspace{-2pt}\mathrm{PAW}}=
\delta_{l_1l_2}\delta_{m_1m_2}\braket{R^a_{j_1}}{R^a_{j_2}}_{\hspace{-2pt}\mathrm{PAW}}\]
i.e. the inner product of partial waves \(\phi^a_i\) is restricted within PAW spheres.parbigskip Expanding the dot products in the penalty functional, we get
\[\begin{split}E_\mathrm{cDFT}=\Lambda&\left(\vphantom{\frac{1}{1}}
\left({1-\hat{u}_{a,x}^2}\right)m_{a,x}^2 + \left({1-\hat{u}_{a,y}^2}\right)m_{a,y}^2 + \left({1-\hat{u}_{a,z}^2}\right)m_{a,z}^2
\right. \\
&\left.\vphantom{\frac{1}{1}}-2\left({\hat{u}_{a,x}\hat{u}_{a,y} m_{a,x}m_{a,y}
+\hat{u}_{a,x}\hat{u}_{a,z} m_{a,x}m_{a,z}
+\hat{u}_{a,y}\hat{u}_{a,z} m_{a,y}m_{a,z}}\right) \right),\end{split}\]
and since we have
\[\frac{\partial m_{a,x}}{\partial D^a_{x,i_1i_2}}=\mu_\mathrm{B}N_{i_1i_2},\qquad
\frac{\partial m_{a,y}}{\partial D^a_{y,i_1i_2}}=\mu_\mathrm{B}N_{i_1i_2},\qquad
\frac{\partial m_{a,z}}{\partial D^a_{z,i_1i_2}}=\mu_\mathrm{B}N_{i_1i_2},\]
we can use the chain rule to get the additions to the atomic Hamiltonians
\[\begin{split}H^a_{x,i_1i_2}=\frac{\partial E_\mathrm{cDFT}}{\partial D^a_{x,i_1i_2}}=2\Lambda\mu_\mathrm{B}&
\left[{\vphantom{\frac{1}{1}}
\left({1-\hat{u}_{a,x}^2}\right)m_{a,x}-\hat{u}_{a,x}\left({\hat{u}_{a,y}m_{a,y}+\hat{u}_{a,z}m_{a,z}}\right)}\right]N_{i_1i_2}, \\
H^a_{y,i_1i_2}=\frac{\partial E_\mathrm{cDFT}}{\partial D^a_{y,i_1i_2}}=2\Lambda\mu_\mathrm{B}&\left[{\vphantom{\frac{1}{1}}
\left({1-\hat{u}_{a,y}^2}\right)m_{a,y}-\hat{u}_{a,y}\left({\hat{u}_{a,x}m_{a,x}+\hat{u}_{a,z}m_{a,z}}\right)}\right]N_{i_1i_2}, \\
H^a_{z,i_1i_2}=\frac{\partial E_\mathrm{cDFT}}{\partial D^a_{z,i_1i_2}}=2\Lambda\mu_\mathrm{B}&\left[{\vphantom{\frac{1}{1}}
\left({1-\hat{u}_{a,z}^2}\right)m_{a,z}-\hat{u}_{a,z}\left({\hat{u}_{a,x}m_{a,x}+\hat{u}_{a,y}m_{a,y}}\right)}\right]N_{i_1i_2},\end{split}\]
- class gpaw.new.constraints.SpinDirectionConstraint(constraint, penalty=0.8)[source]
Spin-direction constraint.
See gpaw/test/noncollinear/test_spin_dir_constraint.py
for how to use the
SpinDirectionConstraint extension.
2D example
https://journals.aps.org/prb/abstract/10.1103/PhysRevB.62.11556
from gpaw.new.ase_interface import GPAW
from ase import Atoms
import numpy as np
a = 6.339
d = 1.331
atoms = Atoms('V3Cl6',
cell=[a, a, 1, 90, 90, 60],
pbc=[1, 1, 0],
scaled_positions=[
[0, 0, 0],
[1 / 3, 1 / 3, 0],
[2 / 3, 2 / 3, 0],
[0, 2 / 3, d],
[0, 1 / 3, -d],
[1 / 3, 0, d],
[1 / 3, 2 / 3, -d],
[2 / 3, 1 / 3, d],
[2 / 3, 0, -d]])
atoms.center(axis=2, vacuum=5)
m = 3.0
magmoms = np.zeros((9, 3))
magmoms[0] = [m, 0, 0]
magmoms[1] = [-m / 2, m * 3**0.5 / 2, 0]
magmoms[2] = [-m / 2, -m * 3**0.5 / 2, 0]
atoms.calc = GPAW(mode={'name': 'pw',
'ecut': 400},
magmoms=magmoms,
symmetry='off',
kpts={'size': (2, 2, 1), 'gamma': True},
parallel={'domain': 1, 'band': 1},
txt='VCl2_gs.txt')
atoms.get_potential_energy()
atoms.calc.write('VCl2_gs.gpw')
# web-page: mag1d.png, mag2d.png
from gpaw.new.ase_interface import GPAW
import matplotlib.pyplot as plt
import numpy as np
calc = GPAW('VCl2_gs.gpw')
dens = calc.dft.densities()
grid_spacing = calc.atoms.cell[2, 2] / 200
nt = dens.pseudo_densities(grid_spacing)
n = dens.all_electron_densities(grid_spacing=grid_spacing)
i = nt.desc.size[2] // 2
x, y = n.desc.xyz()[:, :, i, :2].transpose((2, 0, 1))
uv = n.data[1:3, :, :, i]
m = (uv**2).sum(0)**0.5
u, v = uv / m
fig, ax = plt.subplots()
ct = ax.contourf(x, y, m)
cbar = fig.colorbar(ct)
cbar.ax.set_ylabel('magnetization [Å$^{-3}$]')
ax.quiver(*(a[::3, ::3] for a in [x, y, u, v]))
ax.axis('equal')
ax.set_xlabel('x [Å]')
ax.set_ylabel('y [Å]')
fig.savefig('mag2d.png')
fig, ax = plt.subplots()
x, y = n.xy(1, ..., 0, i)
x, yt = nt.xy(1, ..., 0, i)
j = len(x) // 2
L = calc.atoms.cell[0, 0]
x = np.concatenate((x[j:] - L, x[:j]))
y = np.concatenate((y[j:], y[:j]))
yt = np.concatenate((yt[j:], yt[:j]))
ax.plot(x, y, label='all-electron')
ax.plot(x, yt, label='pseudo')
ax.legend()
ax.set_xlabel('x [Å]')
ax.set_ylabel('magnetization [Å$^{-3}$]')
fig.savefig('mag1d.png')
Experiential:
Theoretical:
DFT: