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 = 4for d-electron systems andLMAXMIX = 6for 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.
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.
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
CartesianKPOINTS 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]. AReciprocal(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 smallCartesiancoordinates so the point truly lies on the Cartesian axis. - Place the plane correctly in Cartesian coordinates. In a
Cartesianlist 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 = -1so 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