Using Phonopy with Wien2k

In this post I will describe how to calculate various phonon properties from Wien2k calculations, using Atsushi Togo's Phonopy program. As an example, I will use lithium hydride (LiH) which has a simple NaCl structure with a lattice constant of 4.01A (optimized in GGA approximation).

1. LiH has an fcc structure, but phonopy requires the P1 symmetry for Wien2k. So we convert LiH.struct to P1 by running x supercell and creating 1x1x1 cell with a type P. This cell should have 8 atoms - 4 Li and 4 H.

2. Create INPHON file with the following content:
ATOM_NAME = Li H
NDIM = 1 1 1
LSUPER = .TRUE.

The number of dimensions should be larger (e.g. NDIM 2 2 2), but for the testing purposes I've set it to 1 1 1. Then run:
phonopy --wien2k=LiH.struct

This creates the DISP file and several struct files with the supercells (in this particular case - LiH.structS, LiH.structS-001, LiH.structS-002. Having a look at the DISP file:
1 1 0 0
5 1 0 0

And having a look at LiH.structS-* files, one can see, that in these files atom 1 (Li) and atom 5 (H) are displaced by 0.01A in x direction. You can change the displacement manually, if you want.

3. Next step is to run two Wien2k calculations for LiH.structS-001 and LiH.structS-002 structure files, respectively, and save the resulting *scf files.

Important! The init for this calculations should begin from symmetry, e.g.:
init_lapw -b -s symmetry -vxc 13 -numk 100 -sp

This is important to start from symmetry because sgroup could reduce the symmetry, which will make phonopy unable to create a proper FORCES file. Also, the RMT spheres can overlap (check it with x nn!) and should be reduced if needed.

4. Once you have LiH001.scf and LiH002.scf from the above calculations, copy them in the directory with INPHON and DISP files, and run the following to generate FORCES:
phonopy --wien2k=LiH.struct -f LiH001.scf LiH002.scf

The order of scf files in this command should match the order of the atoms in DISP! If everything is ok, than the FORCES file will be saved and you can calculate various properties.

5. Backup INPHON to INPHON.1 and modify it for your needs as described in the phonopy manual. It is important that you comment or remove LSUPER=TRUE line. To calculate phonon DOS the new INPHON may look as follows:
ATOM_NAME = Li H
NDIM = 1 1 1
MP = 21 21 21

The phonon DOS can be calculated by invoking the command:
phonopy --wien2k=LiH.struct -p

References
1. Phonopy web page, where you may find the documentation.
2. Using phonopy with VASP (in Chinese, but readable with Google translate for non-speakers).

Parallel compilation of Wien2k on Intel Itanium architecture

Here I briefly describe how to compile Wien2k on the supercomputer with Intel Itanium processors. I suppose that you have read the Wien2k UserGuide and know how to compile the serial version of the program.

What do we need:
Intel Fortran Compiler (ifort), I use version 11.0
Gcc or Intel C/C++ compiler (icc)
Intel CMKL libraries (I have used 10.1.1.019 version)
FFTW 2.1.5 (not 3.0!!) compiled with mpi support
MPI (I use MPICH2)

Compiling serial version
For serial version the compilation options given to siteconfig are the following:

System type: K Linux (Intel ifort 11.x compiler + mkl )
Fortran compiler: ifort
C compiler: cc
Compiler options, BLAS and LAPACK:
O Compiler options: -FR -mp1 -w -prec_div -pc80 -pad -align -DINTEL_VML -traceback
L Linker Flags: $(FOPT) -L/opt/intel/cmkl/10.1.1.019/lib/64 -lguide -pthread
P Preprocessor flags '-DParallel'
R R_LIB (LAPACK+BLAS): $(FOPT) -L/opt/intel/mkl/10.1.1.019/lib/64 -lmkl_lapack -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -lguide -pthread

Note: replace the path to your own one, and mind that lib/64 libraries are needed for Itanium. If you have Intel Xeon or similar processor, use lib/emt64 libraries instead.

Compiling FFTW
This is needed only for the parallel compilation. Download FFTW 2.1.5, unpack it and use the following commands for the compilation:

export CC=icc
export F77=ifort
./configure --enable-mpi --enable-threads
make

Copy the libraries from fftw-2.1.5/fftw/.libs and mpi/.libs somewhere to /opt/fftw or leave them there (as I did).

Compiling in parallel
I have mpich2 installed in /opt/mpich2, so I add /opt/mpich2/bin into the path:

export PATH="$PATH:/opt/mpich2/bin"

Check if mpif90 command works:

$mpif90 -v
mpif90 for 1.0.8p1
Version 11.0

The compilation options will be the following:
Compiler: mpif90

RP RP_LIB(SCALAPACK+PBLAS): -L/opt/intel/mkl/10.1.1.0.19/lib/64 -lmkl_intel_lp64 -lmkl_scalapack_lp64 -lmkl_blacs_lp64 -lmkl_sequential -lmkl -L/opt/mpich2/lib -lmpich -L/home/hartree-fock/fftw-2.1.5/mpi/.libs -lfftw_mpi -L/home/hartree-fock/fftw-2.1.5/fftw/.libs -lfftw
FP FPOPT(par.comp.options): $(FOPT) -I/opt/mpich2/include
MP MPIRUN commando : mpirun -np _NP_ -machinefile _HOSTS_ _EXEC_

That's it! Don't forget to add the necessary lines to .bashrc (LD_LIBRARY_PATH etc).

Magic command to prevent Wien2k crash

Sometimes Wien2k program crashes without any reasons and *.error files are empty or just contain the line like "Error in LAPW0". Then run the following command before running Wien2k:

ulimit -s unlimited

It sets the maximum stacks size for the current user unlimited. For clusters and supercomputers you need to add this line to the script which runs Wien2k.

Wien2k Fedora 10 issue

In Fedora 10 the Wien2k program doesn't work at all because of csh/tcsh compatibility issues, giving the error message like this:

argv1: Subscript out of range

The solution to the problem is described on Wien2k mailing list. You either have to install tcsh package from Fedora 8 (you may find one on rpmfind.net), or do the following:
1. Edit x_lapw script and find the following line (should be around the lines 238—248):

echo "`date`> ($name) $command $argv1[2-]" >> $log

2. Comment it and replace with this one:

echo "`date`> ($name) $argv1" >> $log

Now it should work.

Calculate DOS after band structure calculations in Wien2k

Density of States and band structure are calculated on different k point meshes in Wien2k. One mights notice that you cannot calculate density of states after band structure because you get the following error:

forrtl: severe (24): end-of-file during read, unit 29
...
Image PC Routine Line Source
lapw2 0820F1E4 Unknown Unknown Unknown
lapw2 0820D97D Unknown Unknown Unknown
lapw2 081E1CFB Unknown Unknown Unknown
lapw2 081B840A Unknown Unknown Unknown
lapw2 081CC100 Unknown Unknown Unknown
lapw2 080624C0 fermi_ 67 fermi_tmp_.F
lapw2 080881DD MAIN__ 258 lapw2_tmp_.F
lapw2 0804B906 Unknown Unknown Unknown
libc.so.6 00BA1390 Unknown Unknown Unknown
lapw2 0804B811 Unknown Unknown Unknown
0.215u 0.000s 0:00.21 100.0% 0+0k 0+120io 0pf+0w
error: command /home/hartree-fock/wien2k/source/lapw2 uplapw2.def failed

But there is the way to fix this by running lapw1 without a -band switch. For spin-polarized calculation run:

x lapw1 -up
x lapw1 -dn

If you are running non-spinpolarized case than just run "x lapw1". Also, don't forget about the "-orb" switch if you're using some orbital dependent potentials.

After that you'll be able to calculate DOS without any problems.