WaveTrans: Real-space wavefunctions from VASP WAVECAR file

R. M. Feenstra and M. Widom
Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213

Introduction

This program reads in the binary file WAVECAR that is produced by VASP, and it outputs into a text (ASCII) file the corresponding G values (lattice points in reciprocal space) and their associated plane wave coefficients. From these coefficients, a real-space wavefunction (pseudo-wavefunction) can be constructed in the usual way. Due to its binary format and the ordering of its G values, it can be problematic to directly read and translate the WAVECAR file. The WaveTrans program and its associated WAVECARin routine provides a solution to this problem. The format of the WAVECAR file is:
Record-length #spin components RTAG(a value specifying the precision)
#k-points #bands ENCUT(maximum energy for plane waves)
LatVec-A
LatVec-B
LatVec-C
Loop over spin
   Loop over k-points
      #plane-waves k-Vec
      Loop over bands
         band-energy, band-occupation
      End loop over bands
      Loop over plane waves
         Plane wave coefficient
      End loop over plane waves
   End loop over k-points
End loop over spin
A value of RTAG of 45200 specifies that complex*8 binary format is used in WAVECAR for the coefficients, and a value of 45210 specifies complex*16 format; WaveTrans presently works only for the former case (see 'Limitations' section below). The WaveTrans program reads in the contents of WAVECAR and outputs the GCOEFF file with the appropriate values of the wavevectors and the corresponding plane wave coefficients. For specification of the format of the GCOEFF file, see the comments at the top of the source code. Source code for a FORTRAN 77 version is available in WaveTrans.f and WAVECARin.f (the latter includes the routine vcross for computing cross-products for 3-component vectors). Source code for a FORTRAN 90 version with dynamic allocation of arrays is in WaveTrans.f90. An executable version of the FORTRAN 77 version (made using Gnu gfortran) that should run on any Windows PC is available in the file WaveTrans.exe. For other platforms, simply compile and link the source code on the specific platform. It should be noted that details regarding binary input vary from compiler to compiler. For the case of the ifort compiler (on Linux) the compiler option "-assume byterecl" must be specified so that the record length indicator in the program is interpreted in bytes. Similar types of option(s) may be required for other compilers.

Usage

The program reads in the binary file WAVECAR, that is output from VASP. For the FORTRAN 77 version (WaveTrans.f), the dimension of the various arrays in the program must be sufficiently large; if they are too small then an error message results and the user must then increase the array size (using the PARAMETER statements in the source code) and recompile the program. For the FORTRAN 90 version (WaveTrans.f90), the arrays sizes are set automatically.

Example

From the plane wave coefficients in the file GCOEFF that is output from WaveTran, real-space wavefunctions can be easily constructed (as described in Supplemental material for Phys. Rev. B 87, 041406 (2012)). An example of this procedure is provided by the program WaveFcnPlot, with source code WaveFcnPlot.f and executable code WaveFcnPlot.exe. This sample program produces a single wavefunction, as a function of z, at specified values of x, y, spin, wavevector, and band. The wavefunction is output to WAVEFCN.txt. An example of the output from the program is shown in the posted WAVEFCN.txt file, for a WAVECAR file corresponding to a layer of graphene surrounded by vacuum. The corresponding GCOEFF file is not posted, due to its relatively large size (since it is in text rather than binary format), but it can be obtained simply by running WaveTrans.

A much more efficient way to construct the WAVEFCN.txt file is to write a program that combines the functions of WaveTrans and WaveFcnPlot, without generating the intermediate GCOEFF file. Such a program is also posted, with source code WaveTransPlot.f and executable code WaveTransPlot.exe. This program reads in the WAVECAR file and produces an output file WAVEFCN.txt that is identical to the one produced by the combination of WaveTrans and WaveFcnPlot, but does so using only a small fraction of the CPU time and disk space that the combination requires. The FORTRAN 90 version, WaveTransPlot.f90, offers command line options, e.g. the command "WaveTransPlot -h" produces:

syntax: WaveTransPlot -f file -s spin -k k-point -b band -x coord -y coord
defaults: -f WAVECAR -s 1 -k 1 -b 1 -x 0.0 -y 0.0
inputs: x and y are direct coordinates on axes a1 and a2
output: wavefunction psi(x,y,z) with z direct coordinate on a3 axis

WAVECAR files for spinors (noncollinear magnetism)

All of the above discussion applies to WAVECAR file for collinear magnetism, that is, with either 1 or 2 spin components in the file. VASP also can be used for problems in noncollinear magnetism (e.g. with spin-orbit interaction), in which case the WAVECAR file describes spinor wavefunctions and has a different structure. A FORTRAN 90 program WaveTransSpinor.f90 is available for handling such cases. It is the analog of the WaveTransPlot.f90 program as described above and is used in a similar manner, except that no -s option need be specified in the command string.

Limitations

WaveTrans presently works only for an RTAG value of 45200, corresponding to single precision for the plane wave coefficients. If the WAVECAR file has some other RTAG value then an error message results. To obtain a version that works in double precision, simply change the declaration for the 'coeff' array to be complex*16 rather than complex*8, change the statement that checks the RTAG value, and recompile the program. WaveTrans assumes a full WAVECAR file, as opposed to the reduced-size version created by a Gamma-only calculation.

For interpreting the WAVECAR files, WaveTrans makes an estimate of the maximum magnitude of each component of the G vectors that need to be considered. These estimates may be in error if the particular wavevector that the states are being evaluated at lies outside of the 1st Brillouin zone. Hence, wavevectors should be restricted to the 1st Brillouin zone. (Alternatively, for a wavevector well outside the 1st Brillouin zone, the user can manually increase the values of nb1max, nb2max, and nb3max, as well as npmax, in the program).

Error Messages

For the f77 version, if the array sizes declared in the program are not large enough, then an error message is produced. The user must modify the array sizes in the PARAMETER statement at the top of the program, and re-compile the program.

For both the f77 and f90 versions, the program makes an estimate of the maximum magnitude of each component of the G vectors that need to be considered. For the f90 version an estimate of the maximum number of plane waves is also made, as needed for the allocation of array space. If either of these estimates is incorrect, so that the number of plane waves computed by WaveTrans deviates from the actual number in the WAVECAR file, then an error message results. This was an issue for version 2.0 of WaveTrans.f90 and WaveTransPlot.f90 and for versions 1.1 and previous of WAVECARin.f, in which some of these estimates failed for unit cells that deviated significantly from cubic or orthorhombic. However, these estimates were revised in more recent versions, and they presently seem to be reliable for all cases.

Another situation that can lead to a discrepency between the number of plane waves computed by WaveTrans and that contained in the WAVECAR file has to do with the numerical value of 2m/hbar^2 employed in WaveTrans. This value must be precisely equal to that used in VASP, such the maximum energy cut has identically the same effect in the two programs. The value of 2m/hbar^2 in WaveTrans is presently the same as that in WAVECAR to 9 significant figures. However, a slight adjustment to the final decimal place (or one beyond that) might still be needed. If a user encounters a situation with an error message saying that "the computed number of plane waves is not equal to the input number", then a slight adjustment of the 2m/hbar^2 value (parameter "c" in the DATA statement at the top of the program/routine) is likely needed. If such an adjustment is successfully determined, then please contact feenstra@cmu.edu with the new value, so that it can be incorporated into a future version of WaveTrans.

Cautionary Note

Users should be aware that VASP wavefunctions are based on the projector augmented-wave method, i.e. near the cores of the atoms they differ greatly from true real-space wavefunctions. Thus they are pseudo-wavefunctions, not true wavefunctions. Also, the VASP wavefunctions are not perfectly normalized. Rather, they are orthonormalized with respect to an overlap operator as discussed, e.g., in G. Kresse and D. Joubert "From ultrasoft pseudopotentials to the projector augmented-wave method", Phys. Rev. B 59, 1758 (1999).