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 vector
      Loop over bands
         band energy, band occupation
      End loop over bands
      Loop over bands
         Loop over plane waves
            Plane-wave coefficient
         End loop over plane waves
      End loop over bands
   End loop over k-points
End loop over spin
On the first line of the WAVECAR file, a value of RTAG of 45200 specifies that complex*8 binary format is used for the coefficients, whereas a value of 45210 specifies complex*16 format. WaveTrans presently works only for the former case (see Limitations section below). Note that in the WAVECAR file there is no explicit specification of the G value associated with a given plane-wave coefficient. Rather, the G values must be deduced by the ordering of the coefficients; it is this operation that WaveTrans accomplishes.

The WaveTrans program reads in the contents of the WAVECAR file, and outputs a GCOEFF.txt file with the appropriate G values and corresponding plane-wave coefficients. For full specification of the format of the GCOEFF.txt file, see the comments at the top of the source code. In brief, the GCOEFF.txt file contains, for each wavevector k, band ν, and spin channel, a set of G values and corresponding complex plane-wave coefficients Cν,G(k), from which the real-space wavefunction can be easily constructed (as described in Supplemental material for Phys. Rev. B 87, 041406 (2013)). Within the GCOEFF.txt file, each G value is specified by providing three integers (n1,n2,n3), from which the G value is obtained according to G = n1b1+n2b2+n3b3, where bi are the reciprocal lattice vectors. The set of G values that are contained explicitly in the GCOEFF.txt file (and implicitly in the WAVECAR file) is determined by the maximum plane-wave energy used for the VASP run, i.e. it is all G values for which ℏ2|k+G|2/2m is less than the specified maximum energy. When using WaveTrans, it is best to restrict the wavevector k to the 1st Brillouin zone, as further discussed in the Limitations section below.

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.txt that is output from WaveTran, real-space wavefunctions can be easily constructed (as described in Supplemental material for Phys. Rev. B 87, 041406 (2013)). 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.txt 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.txt 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 of WaveTransPlot, 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
The program outputs a GCOEFF.txt file containing only the coefficients for the selected k-point, band, and spin channel. The z, Real(psi(z)), and Imag(psi(z)) values for the corresponding wavefunction are output to logical unit number 6.

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. Output to the GCOEFF.txt file now contains two sets of wavefunction values for the two spinor components (first one complete set for all G values and then the other), and the output to logical unit number 6 includes two complex wavefunction values for each z value.

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.

Importantly, 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.

In general, WaveTrans computes the number of plane wave coefficients that it is expecting to be supplied with (at each k value), and it compares that with the actual number of coefficients supplied. If these values are not equal, then an error message saying that "the computed no. plane waves not equal to input no." is output, and the program stops. If the "computed no. plane waves" as output by WaveTrans is approximately twice as large as the "input no. plane waves", then probably a Gamma-only computation is being done (see Limitations section above). WaveTrans cannot handle this case.

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/ℏ2 employed in WaveTrans. This value must be precisely the same as that used in VASP, such the maximum energy cutoff has identically the same effect in the two programs. The value of 2m/ℏ2 in WaveTrans is presently equal to 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", with the "computed number" being very nearly equal to the "input number", then a slight adjustment of the 2m/ℏ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).