[−][src]Crate wavecar_rs
Overview
Here is the description of the WAVECAR file structure.
Recordlength #spin components RTAG(a value specifying the precision)
#kpoints #bands ENCUT(maximum energy for plane waves)
LatVecA
LatVecB
LatVecC
Loop over spin
Loop over kpoints
#plane waves, k vector
Loop over bands
band energy, band occupation
End loop over bands
Loop over bands
Loop over plane waves
Planewave coefficient
End loop over plane waves
End loop over bands
End loop over kpoints
End loop over spin
Meta information
The meta information contains two records (see the record concept in fortran's unformatted io action).
The first record only contains
RECL, NSPIN, RTAG
corresponding to recordlength, number of spin components, precision tag. Though these three values are integer, they are stored in float64.
Here is how we treat RTAG:
let prec_type = match RTAG {
45200 => Complex32,
45210 => Complex64,
53300 => Err("Unsupported WAVECAR format: VASP5 with f32")
53310 => Err("Unsupported WAVECAR format: VASP5 with f64")
_ => Err("Invalid WAVECAR format: Unknown VASP version")
};
The second record contains some more information:
NKPTS, NBANDS, ENCUT, LATT, EFERMI
corresponding to number kpoints, number of bands in each kpoint, energy cutoff, lattice vectors in real space (3x3 matrix), fermi energy (in vasp 5 and higher). NKPTS and NBANDS are casted from float64 to integer when reading. ENCUT, LATT and EFERMI are originally float64 value or matrix.
Here LATT is stored in rowmajor, which means in memory, it should be:
LATT[0][0], LATT[0][1], LATT[0][2],
LATT[1][0], LATT[1][1], LATT[1][2],
LATT[2][0], LATT[2][1], LATT[2][2],
Note: Meta information or header is stored in float64 for all the WAVECARs. The main wavefunction
coefficients are stored in either float32 or float64 determined by RTAG
in meta information.
Body
The body content starts at the third record. In this part, band eigen value, fermi weight and band coefficients are stored.
Here is the structure, from the third record:
for ispin in 0..NSPIN {
for ikpoint in 0..NKPTS {
// One record here 1)
for iband in 0..NBANDS {
// One record here 2)
}
}
}

There are 4 + 3*NBANDS values here:
 NPLWS > number of plane waves (number of coefficients) in this kpoint, need to be casted into integer;
 KVEC > kvector for current kpoint, three float64 values;
 The next 3*NBANDS values makes up a NBANDS * 3 matrix, where first two rows mean the
eigen values for each band and the last row means the fermiweight (aka fermi occupation);
band[0].real, band[0].imaginary, fermi_weight[0]; band[1].real, band[1].imaginary, fermi_weight[1]; band[2].real, band[2].imaginary, fermi_weight[2]; ... ... ...

This part dominates the WAVECAR. There are NBANDS records here. Each records contains NPLWS planewave coefficients, and the type of coefficients is either complexf32 or complex64, determined by RTAG;
coeff[0], coeff[1], coeff[2], ... // NPLWS values in total
To sum up, the number of total records in WAVECAR should be 2 + NSPIN * NKPTS * (1 + NBANDS)
.
However, for the calculations that enables spin orbits coupling correction, NSPIN = 1
, but on
each band, two spinor components (equal length) are stored, where the upper and the lower is
are placed in order.
// on each band, spinor up
coeffs[0], coeffs[1], ...
// spinor down
coeffs[0], coeffs[1], ...
This crate only works on little endian machine. You can patch it if you have some feature requests on big endian machine, that's not difficult (changing all the LittleEndian generic parameter to BigEndian shall do the job).
Implementation of transformation of wavefunction in kspace into realspace
We've got the band coefficients, but how to use it? What should we do if we want to visualize it in real space? The relation between kspace and real space is the Fourier transformation. But in which order is the coefficients placed?
FFT grid generation for standard & SOC system
Wavefunction in kspace is a 3D grid, but we must start with 1D to illustrate how the grid is generated:
Suppose we have a 1D grid whose size is ngrid, the 1D fft grid should be
[0, 1, 2, ... ngrid/2] ++ [(1+ngrid/2ngrid), (ngrid/2ngrid), ... 1]
e.g. when ngrid = 11
ret = [0, 1, 2, 3, 4, 5, 5, 4, 3, 2, 1]
And for 3D grid, three directions does the same job
fx = generate_grid(ngrid[0])
fy = generate_grid(ngrid[1])
fz = generate_grid(ngrid[2])
Then combine them (in Fortran, the inner index is the fastest one)
for ifz in fz {
for ify in fy {
for ifx in fx {
fft_grid += [ifx, ify, ifz]
}}}
Each coordinate [ifx, ify, ifz]
in the above will be named with G
in the following.
Now we've a cube in kspace. However the valid wavefunction grid should be a sphere where the
radius is determined by the formula (G + k)^2 / 2 < ENCUT
, where k
is the kvector of current
kpoint.
We've rubbed the kspace cube into a kspace sphere, and the order is consistent with coefficients in WAVECAR. For now, we cannot put the coefficients in the grid directly, and why and how do we do the job will be illustrated in the next section.
Then the arrangement is done for standard and SOC systems. As for the gamma only system, the arrangement is somewhat more complicate.
FFT grid generation for gamma only system
First we need to perform a standard FFT grid generation, then filter the G points, i.e. cut the
sphere and remove half of it. For gamma half of x
direction:
fft_grid.iter()
.filter([gx, gy, gz] {
(gx > 0) 
((gx == 0) && (gy > 0)) 
((gx == 0) && (gy == 0) && (gz >= 0)
})
And for z
direction:
fft_grid.iter()
.filter([gx, gy, gz] {
(gz > 0) 
((gz == 0) && (gy > 0)) 
((gz == 0) && (gy == 0) && (gx >= 0))
})
Reverse Fourier transformation
In VASP's implementation, the real space grid is at least double the size of kspace grid, which
means [rgx, rgy, rgz] = [kgx * 2, kgy * 2, kgz * 2]
. But we also can specify finer grid to
reach better accuracy. Here we follow the VASP's implementation, and take ngrid
as the user
customized realspace grid.
Still remember that there are negative indices in FFT grids? That negative index means count from tail.
let v = [0, 1, 2, 3, 4, 5, 6]; // len = 7
v[0] == 0;
v[1] == 1;
v[1] == 6;
v[2] == 5;
We can transform the negative indices into positive indices via idx_pos = idx_neg.rem_euclid(len)
where len
is length of the dimension where index refers. Go back to our 3D kspace grid, the len
should be the length of corresponding dimension of user customized kgrid.
Now we can put the coefficients on the real gird.
kgrid[0, 0, 0] = coeff[0];
kgrid[1, 0, 0] = coeff[1];
kgrid[2, 0, 0] = coeff[2];
...
If the system is standard of SOC, a simple reverse Fourier transformation would lead to the real space grid.
rgrid = ifft(kgrid)
But for gamma only system, there is still a little job to do:
Because gamma only version only stores half of the sphere, the reverse Fourier transformation
is irfft
(aka complex to real reverse Fourier transformation). This means that we must use a
half of the grid to store the coefficients.
For z
direction, kgrid.shape[2] = rgrid.shape[2]/2 + 1
for fx in gxs {
for fy in gys {
if (fy > 0  (fy == 0 && fx >=0)) continue;
kgrid[fx, fy, 0] = kgrid[fx, fy, 0].conjugate()
}
}
For 'x' direction, kgrid.shape[0] = rgrid,shape[0]/2 + 1
for fy in gys {
for fz in gzs {
if (fy > 0  (fy == 0 && fz >= 0) continue;
kgrid[0, fy, fz] = kgrid[0, fy, fz].conjugate()
}
}
Then rgrid = irfft(kgrid)
can produce the right result.
Visualize the wavefunction in real space
Just take the rgrid
and save in CHGCAR format and visualize it via
VESTA.
Acknowledgement
 Qijing Zheng;
 ExpHP;
 Other guys from the group.
Structs
Wavecar  Main Wavecar struct 
WavecarError  
Wavefunction  Pseudo wavefunction structure. 
Enums
GammaHalfDirection  For gamma only version, there are two gamma half implementations: 
WFPrecisionType  Wavefunction precision type 
WavecarType  WAVECAR type enumeration. 