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 the [math]\displaystyle{ \sigma_x }[/math], [math]\displaystyle{ \sigma_y }[/math], and [math]\displaystyle{ \sigma_z }[/math] spin projections with py4vasp.

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.

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

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 most important arguments of to_quiver are selection, which names the two in-plane spin components and the band index, and supercell, which replicates the plane periodically. Refer to the py4vasp documentation for the remaining arguments.

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

Band-structure calculation using meta-GGA functionals, METAGGA

Files: KPOINTS, CHGCAR, vaspout.h5, PROCAR

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