Computing the spin texture

The spin texture is the momentum-resolved expectation value of the spin, [math]\displaystyle{ \langle\boldsymbol{\sigma}(\mathbf{k})\rangle }[/math], of the Kohn-Sham orbitals. In crystals that lack inversion symmetry, spin-orbit coupling lifts the spin degeneracy of the bands and locks the spin direction to the crystal momentum k. The resulting winding patterns of the in-plane spin, such as the Rashba and Dresselhaus textures, are a direct fingerprint of the underlying symmetry. This page describes how to compute the spin texture as a noncollinear calculation with spin-orbit coupling, sampled on a two-dimensional k plane, and how to read out and plot the [math]\displaystyle{ \sigma_x }[/math], [math]\displaystyle{ \sigma_y }[/math], and [math]\displaystyle{ \sigma_z }[/math] spin projections with py4vasp.

Mind: py4vasp reads the vaspout.h5 file, which VASP writes only when it is compiled with HDF5 support.

The procedure applies to any inversion-asymmetric system once the relevant spin-split bands and the appropriate k plane have been identified.

Step-by-step instructions

Step 1: Compute the self-consistent ground state

Set up a self-consistent noncollinear calculation with spin-orbit coupling and run it with the vasp_ncl executable. The relevant INCAR tags are:

ISMEAR = 0
SIGMA = 0.05
LSORBIT = .TRUE.       ! spin-orbit coupling (also sets LNONCOLLINEAR)
MAGMOM = 9*0           ! three components per atom; here three atoms, no initial moment
SAXIS = 0 0 1          ! spin-quantization axis (default)
GGA_COMPAT = .FALSE.   ! improves the numerical precision of GGA with SOC
LASPH = .TRUE.         ! aspherical one-center contributions
LMAXMIX = 4            ! d electrons: write the charge density to higher l
LCHARG = .TRUE.        ! write the converged density for the restart

LSORBIT=.TRUE. switches on spin-orbit coupling and automatically sets LNONCOLLINEAR=.TRUE., so the calculation treats the full 2×2 spin density. In the noncollinear case, MAGMOM takes three components per atom; the example uses a nonmagnetic starting guess. The spin-quantization axis is fixed by SAXIS, and all spin projections are reported in this basis. Keeping symmetry switched on for the ground state speeds up the self-consistent run.

Important: Set LMAXMIX = 4 for d-electron systems and LMAXMIX = 6 for f-electron systems so that the one-center PAW charge densities are written to the CHGCAR file for the non-self-consistent restart.

Step 2: Compute the band structure and identify the Rashba-split bands

Before sampling a plane, compute the one-dimensional band structure along a high-symmetry path and locate the spin splitting. Follow the procedure in Band-structure calculation using density-functional theory to set up and plot the band structure; restart from the converged density of Step 1.

Choose a path that passes through the time-reversal-invariant point where you expect the splitting. The Rashba signature in the band structure is:

• The spin-split bands are Kramers-degenerate at the time-reversal-invariant point and split as k moves away from it.
• The dispersion shows the characteristic Rashba shape: the band extremum is displaced off the high-symmetry point, so the band develops a local extremum at the point and a ring of true extrema around it (the "camelback" for a conduction band).

This identification tells you which bands carry the texture and on which point to center the planar mesh in the following steps. The spin texture is a symmetry effect and is well-defined even for small splittings, so a small splitting does not prevent resolving the texture. However, if the Rashba splitting is smaller than expected, it may indicate a problem with the setup, such as an inaccurate geometry or a too-small bandgap.

Example band structure along a high-symmetry path through the time-reversal-invariant point. The spin-split pairs are degenerate at the point and split away from it, with the band extrema displaced off the point. This example is bulk BiTeI along H-A-L.

Step 3: Choose the k plane from the crystal symmetry

The spin texture is a two-dimensional quantity, so the non-self-consistent run samples a planar k mesh rather than a path. Use the band structure to identify the band extremum and to confirm that a Rashba splitting is present there. It does not necessarily have to be the largest splitting in the band structure.

The planar mesh has a single point along the direction normal to the plane. Depending on the symmetry of the material, the mesh therefore takes the form n n 1, n 1 n, or 1 n n, with the shift selecting the value of the fixed component.

To find the relevant plane for a given material:

• Identify the symmetry responsible for the texture. Rashba and Dresselhaus splittings both require broken inversion symmetry; the spin-momentum locking is organized around a high-symmetry point or axis.
• Orient the plane perpendicular to the relevant symmetry or polar axis. For a polar crystal with a unique polar axis, the Rashba spin-orbit field lies in the plane perpendicular to that axis, so sample a k plane perpendicular to it.
• Center the plane on the band extremum identified in Step 2, and set the KPOINTS shift so that the fixed component matches that extremum.
• Plot the two in-plane components perpendicular to the plane normal. A nonzero in-plane winding (a vortex or anti-vortex) is the Rashba signature, and the out-of-plane component is close to zero near the center.

Step 4: Run the non-self-consistent calculation on the planar mesh

Create a directory for the planar run and copy the input files together with the converged density:

mkdir -p texture cp INCAR POSCAR POTCAR CHGCAR texture/. cd texture

Provide the planar k mesh in the KPOINTS file. The shift selects the plane that contains the extremum identified in Step 2. If the band extremum is at the [math]\displaystyle{ \Gamma }[/math] point, no shift is needed. Otherwise, choose the shift so that a mesh point coincides with the extremum: shift the out-of-plane component to select the plane that contains the extremum, and shift the in-plane components if the extremum is not at the projected zone center.

The following is an example for bulk BiTeI, where the extremum (the A point) lies at the zone boundary along the out-of-plane direction, so the third component is shifted to 1/2:

Planar mesh (example: BiTeI A-plane)
0
Gamma
11 11 1
0 0 0.5

Add the following tags to the INCAR file from Step 1 for a restart at fixed density:

ICHARG = 11        ! restart from CHGCAR, density held fixed
ISYM = -1          ! keep the full planar mesh
LORBIT = 11        ! write the spin projections
LWAVE = .FALSE.
LCHARG = .FALSE.

ICHARG = 11 reads the converged charge density and keeps it fixed. ISYM = -1 switches off symmetry, which is essential here: with symmetry on, VASP reduces the planar mesh to the irreducible wedge, but the quiver plot requires the full planar mesh. LORBIT = 11 writes the site- and orbital-resolved spin projections that py4vasp reads.

To refine the texture, increase the mesh density (for example to 21×21×1); the plane always spans the full Brillouin zone, and the shift only selects the fixed component.

Step 5: Plot the spin texture with py4vasp

py4vasp draws the in-plane spin as a quiver plot. Run the following in a Python notebook:

import py4vasp
calc = py4vasp.Calculation.from_path("/path/to/calculation")
conduction_band = 29  # replace with the index of the conduction band in your calculation
calc.band.to_quiver(f"sigma_x~sigma_y(band={conduction_band})")

The arguments of to_quiver are selection, which names the two in-plane spin components and the band index, supercell, which replicates the plane periodically, and normal which rotates the cell in the plan. Refer to the py4vasp documentation for the details.

Example of the to_quiver output over the full Brillouin zone: the in-plane spin arrows wind around the band extremum, and the threefold pattern reflects the crystal symmetry. This example is bulk BiTeI on a plane through the top face of the Brillouin zone.

Verification

A spin-texture plot is hard to judge in isolation. Probe a few explicit k points near the center of the plane and read out the spin components directly to isolate the physics from any meshing or plotting artifact.

Probe explicit k points near the center

Pick a few k points a small distance from the band extremum and read out the spin components directly. First, provide an explicit Cartesian KPOINTS list with small displacements along [math]\displaystyle{ \pm x }[/math] and [math]\displaystyle{ \pm y }[/math]:

Explicit k points near the band extremum
4
Cartesian
  0.01  0.00  0.00  1
 -0.01  0.00  0.00  1
  0.00  0.01  0.00  1
  0.00 -0.01  0.00  1

Set the third component so that the points lie on the plane of the extremum, as in Step 4, and mind the Cartesian-unit caveat noted under Recommendations. Run this with the same restart INCAR as Step 4.

Then read the spin components and print one of them for the band of interest:

import py4vasp
calc = py4vasp.Calculation.from_path("/path/to/calculation")
spin = calc.band.read("x, y, z")  # keys "x","y","z"; each array has shape [k point, band]
# read() returns a 0-based array, while the band index in the to_quiver selection
# string is 1-based, so subtract one here
print(spin["x"][:, conduction_band - 1])  # sigma_x at every k point for the conduction band

The expected pattern depends on the type of spin-orbit coupling:

• Rashba (polar, [math]\displaystyle{ C_{nv} }[/math]): the spin is tangential, perpendicular to k in the plane, winding in a single circular sense, with [math]\displaystyle{ \sigma_z\approx0 }[/math] near the center. For a point along the Cartesian x axis only [math]\displaystyle{ \sigma_y }[/math] is nonzero, and along y only [math]\displaystyle{ \sigma_x }[/math].
• Dresselhaus (bulk-inversion asymmetry, for example zinc-blende): the pattern is not a simple circle; the spin direction follows the cubic angular dependence, so the perpendicular-to-k rule does not hold in the same way.

The idea of sampling near the center and checking the component signs against the expected symmetry is general; only the expected pattern changes with the type of coupling.

Recommendations and advice

• Use the noncollinear executable. Spin-orbit coupling requires vasp_ncl; the standard or gamma-only executables do not include it.
• Mind the k-point units. A Cartesian KPOINTS list is given in units of [math]\displaystyle{ 2\pi/\text{SCALE} }[/math], where SCALE is the POSCAR scaling factor. A value of 0.05 then corresponds to [math]\displaystyle{ 0.05\cdot2\pi\approx0.31 }[/math] Å^-1, which can fall far outside the linear Rashba region and into a regime of warping with a spurious [math]\displaystyle{ \sigma_z }[/math]. A Reciprocal (fractional) list avoids the [math]\displaystyle{ 2\pi }[/math] factor, but its axes are the reciprocal lattice vectors, which are non-orthogonal for a hexagonal cell. There, [math]\displaystyle{ \mathbf{b}_1 }[/math] is 30° off the Cartesian x axis, so a fractional [math]\displaystyle{ (\delta, 0, *) }[/math] point is not along x and its tangential spin acquires a nonzero [math]\displaystyle{ \sigma_x }[/math] even for an ideal Rashba state. To test the "only [math]\displaystyle{ \sigma_y }[/math] along x" rule, use small Cartesian coordinates so the point truly lies on the Cartesian axis.
• Place the plane correctly in Cartesian coordinates. In a Cartesian list the third coordinate is also scaled by [math]\displaystyle{ 2\pi/\text{SCALE} }[/math]. To sit on a plane at fractional [math]\displaystyle{ k_z=1/2 }[/math] use [math]\displaystyle{ k_z=0.5/c }[/math] rather than 0.5.
• to_quiver needs a full regular mesh. The quiver routine reads the grid dimensions from the mesh metadata and rejects an explicit k list. The plotted plane therefore always spans the whole Brillouin zone; refine the texture by increasing the mesh density, not by narrowing the window. The routine plots the in-plane components only and does not color the arrows by [math]\displaystyle{ \sigma_z }[/math].
• Keep symmetry off for the spin-texture calculation. Set ISYM = -1 so that all k points are kept; the full mesh is needed for plotting.

Related tags and articles

Band-structure calculation using density-functional theory

Files: KPOINTS, vaspout.h5, PROCAR

Tags: LSORBIT, LNONCOLLINEAR, SAXIS, MAGMOM, LORBIT, ISYM, ICHARG, LMAXMIX, GGA_COMPAT, LASPH