I'm Bk Ryu from South Korea.
My major is the theoretical solid state physics.
I'm a specialist in predicting physical properties from computer simulations and thermoelectric things, as well as the electronic structure of defects in solid systems.
Here I'd like to share my tips how to use VASP efficiently.
When you copy my web page content, please be sure to make a link to my site (20170519). [Copyright 2017, 2016, 2015, 2014 by Byungki Ryu]
 VASP manual site
http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html
 LINK for VASP tips (update: 20150427)
David Karhánek's Homepage ( http://homepage.univie.ac.at/david.karhanek/downloads.html )
Nano Group from Budapest ( http://wiki.kfki.hu/nano/VASP_relevant_resources )
WaveTrans: Realspace wavefunctions from VASP WAVECAR file (Dr. R. M. Feenstra and Dr. M. Widom) ( http://www.andrew.cmu.edu/user/feenstra/wavetrans/ )
JR Kitchin, Modeling materials using density functional theory ( http://kitchingroup.cheme.cmu.edu/dftbook/dft.html ) Skelton's tips: Phonons & Phonopy: Pro Tips (2014) ( http://www.slideshare.net/jmskelton/phononsphonopyprotips2014 ) Phonons & Phonopy: Pro Tips (2015) ( http://www.slideshare.net/jmskelton/phononsphonopyprotips2015 ) VASP: Some Accumulated Wisdom ( http://www.slideshare.net/jmskelton/vaspsomeaccumulatedwisdom ) VASP And Wannier90: A Quick Tutorial ( http://www.slideshare.net/jmskelton/vaspandwannier90aquicktutorial ) Wannier90: Band Structures, Tips and Tricks ( http://www.slideshare.net/jmskelton/wannier90bandstructurestipsandtricks )
 Atomic structure optimization : to obtain optimized structures (update: 20130824)
 In many cases, the atomic positions are drifted if one use default option for structure relaxations. In other words, atomic forces are hardly reaching small values due to the uncertainty in calculation.
 To obtain optimized geometry, it is suggested to use the following tag in INCAR.
ADDGRID = T
 Or you can increase calculation precision by lowering "EDIFF" or by increasing "NELMIN"
It will increase the precision of the force and thereby you wil obtain the optimized structure.
NELMIN = 6

Importance of kspace projection (LREAL=FALSE)
 Not same energies for same structures, LREAL=Auto or LREAL=TRUE
 Sometimes, the same structures can give different total energies if either LREAL=TRUE or LREAL=Auto are used.
For example, one makes same structure but the atomic basis are equally displaced by constant vector.
Many cases it is ok, because the energy difference is small, less than few meV/atom.
However, when one calculates the formation energy, the error can be large if the supercell size is very big.
If it happens, you should recalculate the total energy using kspace projection scheme. (LREAL = False)
 When you have to calculate very precise calculation for phonon mode or for small barrier calculations,
it is strongly recommended to use the "LREAL=F" tag, with other high precision settings even though the supercell size is very large.

Unit of densityofstates (DOS) in DOSCAR file (update: 20130824)
 The unit of DOS is the states per electrons. If you run the single interation, you will always obtain "states/eV"
 However, sometimes you will get different unit in DOSCAR. In many cases, it happens becasue you run the molecular dynamics simulation or structure relaxation.
Then the DOSCAR contatins the sum or average value of DOSs. Note that the NSW value is important while the number of interation is not.
 Densityofstates (DOS) on the vacant site (DOSCAR on the vacant site) (from here,
http://cms.mpi.univie.ac.at/vaspforum/forum_viewtopic.php?4.1692 site broken )
 One can project the electronic wavefunction on arbitrary positions. It is very helpful to find the integrated charge for given spheres without postprocessing of charge density files (CHG or CHGCAR)
 To do that, you should add following lines in POSCAR file, with the use of proper RWIGS, NPAR, LORBIT tags in the INCAR file.
Empty (spheres)
1
0.0000000000000000 0.000000000000000 0.0000000000000000

Drawing local density of states (LDOS) (update: 20131027)
 Local density of states is the real space distribution of the wavefunction square, i.e. charge density for given energy range. Following below steps to obtain the LDOS
 Obtain wavefunction > obtain all partial charge density, rho(n,k,x,y,z) > You can obtain energy eigenvalue from EIGENVAL > Now you know rho(E(n,k),x,y,z)
> Now average over rho for y,z. to obtain rho(E,x) > Now draw rho(E,x) with smearing.
 Note that you need very many kpoint. I know and you know that it is the time consuming job.
 Localpotential (LOCPOT)
 Frequently, one should check the difference between LOCPOT files
 In this case, it is a good idea to fix the number of chargedensity or localpotential grids.
 Related tags are NGX, NGY, NGZ. Maybe sometimes with NGXF, NGYF, NGZF
 Phonon calculations
 Phonopy code is the very useful tool for phonon calculations.
 Download from here: http://phonopy.sourceforge.net/
 Good texts to read:
ftp://ftp.heanet.ie/.../introductionphononcalc.pdf (broken site)
http://icms3.weebly.com/uploads/3/5/9/0/3590130/version1.pdf

 Dielectric function for metallic system
 Dielectric functions are composed of interbandtransitionterm (bound electron term) plus intrabandtransitionterm (free electron term: Drude term)
 Although you can obtain inerbandtransition term from VASP, you should calculate intrabandtransitionterm by yourself especially for metallic system.
 It means that you need plasma frequency from the band structure.
There are two possible method to calculate the plasma freqeuncy using vasp code.
 Dielectric constant for dielectrics with clampedion (update: 20160128)
 Using the LEPSILONtag, one can calculate the ionclamped dielectric constant as written in VASP manual [LINK].  Optimizing structure  Calculating wavefunction  Turning on LEPSILONtag (LEPSILON = TRUE) and commenting out NPARtag or no NPAR tag in the INCARfile.
 Dielectric constant for dielectrics with relaxedion (update: 20160128)
 Using the LEPSILONtag & IBRIONtag, one can calculate the ionclamped dielectric constant as written in VASP manual [LINK].  Optimizing structure. Note that ionic forces should be very small because of the force constant calculations.  Calculating wavefunction  LEPSILON = TRUE / no NPAR tag / IBRION = 5,6,7,8 (one of them)
 Prediction of plasma frequency (update: 20131027)
 The plasma frequency I mention here is not the 'Screened plasma frequency'.
You shoud distinguish between just plasma frequency and screened plasma frequency.
The screened plasma frequency is the frequency where the real part of the dielctric function (interband term + intraband term) becomes zero.
And the screened plasma frequency is smaller than the plasma frequency.
 Using BoltzTraP code, you can obtain electrical conductivity over electon scattering time (SIGMA/TAU) within constant relaxation time approximation.
Then, SIGMA/TAU can be transformed to plasma frequency.
You shoud be careful when calculating plasma frequency.
As the plasma frequency is obtaiend from the intraband transition, I mean the group velocity of the electron,
the kpoint sampling is very important. You need sufficiently many kpoints to obtain correct band dispersion.
It is recommended to test the energy convergence with increasing number of kpoints.
 You may also obtain the plasma frequency within vasp code. With the input tag of "LOPTICS=TRUE" you can find the "plasma frequency" in the OUTCAR file.

Electrical conductivity using VASP code (update: 20131113, 20150417)
 Within semiclassical theory, one may obtain the electrical conductivity from the band structure by calculating the intraband structure.
However, when you use the semiclassical theory, you should get the relaxation time from outside or you should calculate the relaxation time.
If you already know the electron scattering time or relaxation time, you can get the electrical conductivity using VASP.
 There are two ways to get the electrical conductivity. One is using "BoltzTraP" code.
And the other thing is to read the OUTCAR file from VASP code.
Using LOPTICS = TRUE tag and RTIME = x.xx (femtotsec) in INCAR with sufficiently many kpts, you can finally found the reliable value of electrical conductivity tensor in OUTCAR from VASP.
Unit of electrical conductivity per relaxation time is 10^22 [S/m/sec].
For example, if the tensor is [ [0.020 0.000 0.000] [0.000, 0.020, 0.000] [0.000, 0.000, 0.020]], then Sigma/Tau = 0.020 x 10^22 S/m/sec = 2 x 10^20 S/m/sec
For recent VASP version of 5.3.3, you can set the relaxation time ( RTIME in unit of femtosec ) and obtain the electrical conductivity tensor in unit of mega S/m written in OUTCAR. // 20150417
And the Energy dependent conductivity is also written in the "vasprun.xml" file. // 20150417
 MAGMOM tag in INCAR
 You may reduce the sentence length of MAGMOM tag as shown below,
MAGMOM = 1 1 1 1 1 ==> MAGMOM = 5*1
MAGMOM = 1 1 3 3 4 ==> MAGMOM = 2*1 2*3 4
It really saves your time to make INCAR file for spin poloarized calculation of large supercell
 Many body perturbation theory (MBPT) : GW and BSE in VASP (I tested MBPT calculations with vasp version of 5.2.xx)
 MBPT calculations are very time consuming and very memory consuming.
 BSE tips in other site : http://cms.mpi.univie.ac.at/vaspforum/forum_viewtopic.php?4.12605 (site broken)
 Calculation flow:
(1) DFT(Exc=LDA or GGA) calculations (or DFT+U, HybridDFT) for wavefunctions (WAVECAR) and their derivatives (WAVEDER). You need sufficient many unoccupied bands.
(2) GW calculations (ALGO=GW0, scGW0, scGW)
(3) BSE calcualtoin (ALGO=BSE) : although you can not find the BSE tag in your manual, it works. And check the "vasprun.xml" which contains frequency dependent dielectric functions.
 Related tags:
LOPTICS = T for derivative of wavefunction (see manual)
NBANDS number of total bands including occupied and unoccupied, convergence test is needed (see manual)
ENCUT convergence test is needed (see manual)
CSHIFT smearing value for optical response results (freq. dependent dielectric function) (see manual)
ALGO (=GW0, scGW0, scGW, BSE) one shot G0W0 tag (GW0), W fixed GW tag (scGW0), fully self consistent GW (scGW), BSE calculation tag (see manual)
NOMEGA (see manual)
NBANDSO number of occupied bands you considered for BSE calculations
NBANDSV number of unoccipied (virtual) bands you considered for BSE calculations
 Charge analysis (update: 20150420)
 Bader Charge
You may analyze the atomic charge density distribution by using Bader Analysis. [http://theory.cm.utexas.edu/henkelman/code/bader/]
The concept of Bader Charge is based on the space separation based on the minimum charge density positions.
 CHGCAR or CHG file
You can directly analyze the atomic charge density using CHGCAR file.
Or as mentioned in above (#6), you can get the atomic charge using DOSCAR or OUTCAR.
 Unit conversion for VASP phonon frequency from density functional perturbation theory (DFPT) results, Unit of phonopy (update: 20150624)
 [VASP] When one run DFPT calculations or BSE calculations, one obtain the vibration frequency in OUTCAR or vasprun.xml files.VASP.
freq (THz) : ang. moment w (2pi THz) : 1/wavelength (1/cm) : energy (meV)
= 1 : 2 pi : 10^10/c(speed of light: c) : 10^15xh(planck const: h) = 1 : 6.238186 : 33.35641 : 4.135669
 [phonopy] the default unit for frequency is THz. Therefore you can obtain phonon energy by producting planckconstant (h)
1 THz = 4.135669 meV
 example of PbTe phonon (VASP + phonopy, cubic lattice, lattice parameter of a=6.567843.

First principles calculations of Lattice thermal conductivity
 Thermal conductivity (kappa_tot) is a sum of electrical term and lattice term.
 Electrical term (kappa_elec) = (Lorenz Number) x ( Elec. conductivitiy ) x (Temp.), where you can obtain Lorenz number and Elec. conductivity from the BoltzTrrap (electron boltzmann transport equation)
 Lattice therm (kappa_latt) = (kappa_phonon) can be obtained by performing phono3py calculations, which calculates the transition rate of phononphonon scattering. (phono3py LINK)

Using phono3py code with VASP (20151221)
(0) prepare the optimized supercell with the VASP parameter to be used in force calculations.
The force should be less than 1.0e08 eV/angstrom. Of course, you should optimized the structure considering the kpoint mesh you will use in (1).
WAVECAR file is needed to reduce computational resources for step (2).
(1) generate POSCAR files with pairdistance configurations. There might be one to ten thousands POSCAR files.
phono3py d dim="4 3 3" c POSCARunitcell cutoff_pair_distance="4.4"
supercell size, unitcell file, cutoff pairdistance in angstrom
(2) calculate force for pairdistance configurations using VASP.
Make directories for configurations of POSCAR00365.
Copy INCAR, KPOINTS, POTCAR POSCAR, and WAVECAR.
Using WAVECAR preapared from (0) will be useful to reduce the computational cost.
In my case, I used the gamma point only but with the large supercell containing 200400 atoms.
Energy cutoff is also known to be sensitive to the lattice thermal conductivity tensor.
For example, when the monoclinic unitcell volume is optimized within 400 eV and internal coordinate is optimized within 300 eV.
then there could be nonzero offdigonal term in lattice thermal conductivity tensor.
(3) collect vasprun.xml file
phono3py cf3 disp{00001..00999}/vasprun.xml (or) phono3py cf3 disp*/vasprun.xml
disp_fc3.yaml, POSCARunitcell files should be in same directory.
(4) create fc2.hdf and fc3.hdf
phono3py dim="4 3 3" c POSCARunitcell
"This step is not madatory, but you can avoid calculating fc2 and fc3 at every run tim." from "http://phonopy.sourceforge.net/phono3py/workflow.html"
The file size of fc3.hdf is very large. It is about 1 to 10 GB.
(5) calculate thermal conductivity
(5.1) split thermal 3phonon process configurations
phono3py fc3 fc2 dim="4 3 3" v mesh="11 11 11" c POSCARunitcell br thm wgp
qmesh grid, write grid point for split calculation
The computational time for 3phonontransition rates are very time consuming when qmesh grid is fine.
Thermal conductivity, especially for low temperature region, qmesh grid affects the lattice thermal conductivity.
Therefore it is recommended to split the calculation set and calculate the transition rate separately using many cpu cores.
After run the commands in (5.1), you will get "grid_addressm111111" and "ir_grid_points.yaml"
(5.2) calculate 3phonon transition rate for each q grid point
phono3py fc3 fc2 dim="4 3 3" v mesh="11 11 11" c POSCARunitcell br thm write_gamma gp="3"
Check q grid point number (gp) from ir_grid_points.yaml and write correct gp.
The calculation time is dependent on qmesh size.
Each gppoint calculation time is varying from 10 minutes to 10 days depending supercell size and qmesh size.
Note that there are many gp if qmesh size is fine.
(5.3) sum up thermal conductivity
phono3py fc3 fc2 dim="4 3 3" v mesh="11 11 11" c POSCARunitcell br thm read_gamma > z_TC_def.txt
Now you will get the lattice thermal conductivitiy file, "z_TC_def.txt"
Type "tail n 150 z_TC_def.txt" and find the lattice conductivity tensor as a function of temperature.
(*) Useful tag: isotope, mv, bmfp for isotope effect, mass variance effect, boundary scattering effect.
(**) Auxilary tools
kaccum and gaccum commands for accumulated lattice thermal conductivity and 3ph transition rate, respectively.

Band structure using VASP+Wannier90 (20160405) We can draw band structure by calculatinng Wannier Interpolation using Wannier90. You need wannier90 compiled VASP. Although I don't know well about wannier90, I want to share my knowledge to prevent other's failures Simple method (11) Preconverge wavefunction (12) Turn on LWANNIER90=T tag and run wan90 compiled VASP. Then Maximally Localized Wannier Functions (MLWFs) are generated with wannier90.win, wannier90.mmn, wannier90.eig, wannier90.out files. Note that, in my case, the vasp run is failed when I use NPAR tag. If you have same problem, comment out the NPAR tag in the INCAR file. (20170324) (13) ReRun wannier90. wannier90.x wannier90 (14) Interpolate wannier functionals: first modify the wannier90.win file. For details, see the tutorial file. wannier90.x wannier90 GW band structure from wannier90 I will explain how to make GW band structure. You considered diamond Si and there are two Si atoms in primitive cell. Then there are 8 electrons. For GW calculations you may need many number of bands, about 96 or larger than it. Please note that NPAR should be tagged out or equal to number of nodes(cores). And you may know that NBANDS is dependent on NPAR tag. First check NBANDS without NPAR tag. Also you don't need too many wannier functions. It is a good idea to reduce the number of bands. Actually the number of bands for wannier90 is very sensitive to the band structure. So I recommend to check the suitable number of bands by just running simple DFT. Here I assume that the optimal number of wannier bands is 4 per Si atom (8 per primitive cell). In shortly, you need to determine "NBANDS" for VASP and "NUM_WAN" for wannier90. Here I will use NBANDS=96 (I have 12 core nodes, so NBANDS should be divived by 12) Here I will use NUM_WAN = 8 (the number of wannier functions are related to the symmetries) (11) DFT calculations with default NBANDS (12) Optical calculations (LOPTICS=TRUE, NBANDS=96 ) (13) GW calculations with wannier90 (ALGO=GW0, LWANNIER90=TRUE, no NPAR related tag) If you do not set wannier90.win file, then wannier90.win will be generated. But it is too heavy to play with it since there are too many BANDs. Make wannier90.win file (or modify wannier90.win file) as follows. (The default wannier90.win can be obtained by running without wannier90.win file. GW calculation is too time consuming. Therefore just run DFT with LWANNIER90=T tag) num_wann = 8 num_bands = 8 exclude_bands : 996 (14) Run wannier90 After getting wannier90.win, wannier90.eig, wannier90.mmn, rerun the wannier90.x cmd> wannier90.x wannier90 Now you wannierize and obtain MLWFs. (15) Intepolate First modify or generate or make suitable wannier90.win or NAME.win. And run by typing "wannier90.x wannier90" or "wannier90.x NAME". Now you can interpolate for DOS, BANDS, or FermiSurface. IMPORTANT TIP for wannier90 with VASP (20161102): you should comment out the NPAR tag when the calculation is crashed. When I ran the vasp calculation for Mg2Si primitive cell, my calculation is crashed with following error: internal error in GENERATE_KPOINTS_TRANS: G vector not found. By removing NPAR tag, I did succeed to calculate the wannier90 with VASP. Phonon thermal conductivity of low dimensional structure (2016/04/25) Please be careful that if there is vacuum, I mean your structure is 1d or 2d structures with vacuum, you should normalize the thermal conductivity tensors. Also be careful that there will be dipole interactions between adjacent supercell and the force from displacement supercell will be dependent on the lateral size of supercell. For example, I tested 1D carbon (C) atomic chain along z direction and the CC distance is set to 1.28 angstrom. I considered 1x1x10 supercells containing 10 C atoms with various a=b values, 10, 14, 20, and 25 angstroms. I find that the thermal conductivity value of zz component times a square seems to be converged as a becoming larger from 10 to 25 (k_zz x a^2 = 2836, 2968, 3136, 3215 Wxm/K at 300K, 4560, 5181, 5281, 4685 Wxm/K at 500 K. It implies that the cross section area should be checked carefully because the phono3py code assumes that one's structure is bulk. The reason why the k_zz is not exactly same is that there is dipole when single C atom is displaced and it makes large cross sectional interactions between adjacent supercells. I think that if we consider the neutral atomic chain which has negligible dipole interaction between supercell, the k_zz will be converged more fastly.
World wide visitors
