Tips for vasp (tips4vasp)

Original link for google website (click below figure)

Notice

Welcome all on-line and off-line visitors for coming my lab.

Feel free to contact me. I can support you for the travel expenses as well as accommodation as guest house of KERI.

What you should do is to give a seminar or presentation about your work :)

-------------

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 (2017-05-19). [Copyright 2017, 2016, 2015, 2014 by Byungki Ryu]

Tips for vasp users (written by Byungki RYU)

(copied from my original page https://sites.google.com/site/cta4rbk/home/tips4vasp)

1. VASP manual site

   1. http://cms.mpi.univie.ac.at/vasp/vasp/vasp.html

2. LINK for VASP tips (update: 2015-04-27)

   1. David Karhánek's Homepage ( http://homepage.univie.ac.at/david.karhanek/downloads.html )

   2. Nano Group from Budapest ( http://wiki.kfki.hu/nano/VASP_relevant_resources )

   3. WaveTrans: Real-space wavefunctions from VASP WAVECAR file (Dr. R. M. Feenstra and Dr. M. Widom) ( http://www.andrew.cmu.edu/user/feenstra/wavetrans/ )

   4. JR Kitchin, Modeling materials using density functional theory ( http://kitchingroup.cheme.cmu.edu/dft-book/dft.html )

   5. Skelton's tips: Phonons & Phonopy:

      1. Pro Tips (2014) ( http://www.slideshare.net/jmskelton/phonons-phonopy-pro-tips-2014 )

      2. Phonons & Phonopy: Pro Tips (2015) ( http://www.slideshare.net/jmskelton/phonons-phonopy-pro-tips-2015 )

      3. VASP: Some Accumulated Wisdom ( http://www.slideshare.net/jmskelton/vasp-some-accumulated-wisdom )

      4. VASP And Wannier90: A Quick Tutorial ( http://www.slideshare.net/jmskelton/vasp-and-wannier90-a-quick-tutorial )

      5. Wannier90: Band Structures, Tips and Tricks ( http://www.slideshare.net/jmskelton/wannier90-band-structures-tips-and-tricks )

3. Atomic structure optimization : to obtain optimized structures (update: 2013-08-24)

   1. 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.

   2. To obtain optimized geometry, it is suggested to use the following tag in INCAR.

   3. ADDGRID = T

   4. Or you can increase calculation precision by lowering "EDIFF" or by increasing "NELMIN"

   5. It will increase the precision of the force and thereby you wil obtain the optimized structure.

   6. NELMIN = 6

4. Importance of k-space projection (LREAL=FALSE)

   1. Not same energies for same structures, LREAL=Auto or LREAL=TRUE

   2. 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 re-calculate the total energy using k-space projection scheme. (LREAL = False)

   3. 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.

5. Unit of density-of-states (DOS) in DOSCAR file (update: 2013-08-24)

   1. The unit of DOS is the states per electrons. If you run the single interation, you will always obtain "states/eV"

   2. 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.

6. Density-of-states (DOS) on the vacant site (DOSCAR on the vacant site) (from here, http://cms.mpi.univie.ac.at/vasp-forum/forum_viewtopic.php?4.1692 site broken )

   1. 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)

   2. 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

7. Drawing local density of states (LDOS) (update: 2013-10-27)

   1. 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

   2. 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.

   3. Note that you need very many k-point. I know and you know that it is the time consuming job.

8. Local-potential (LOCPOT)

   1. Frequently, one should check the difference between LOCPOT files

   2. In this case, it is a good idea to fix the number of charge-density or local-potential grids.

   3. Related tags are NGX, NGY, NGZ. Maybe sometimes with NGXF, NGYF, NGZF

9. Phonon calculations

   1. Phonopy code is the very useful tool for phonon calculations. Download from here: http://phonopy.sourceforge.net/

   2. Here are some good texts to read:

      1. ftp://ftp.heanet.ie/.../introduction-phonon-calc.pdf (broken site)

      2. http://icms3.weebly.com/uploads/3/5/9/0/3590130/version1.pdf

10. IR-active mode calculations

    1. You can calculate the IR-active mode intensities using the following code.

    2. Download from here: http://homepage.univie.ac.at/david.karhanek/downloads.html

11. Dielectric function for metallic system

    1. Dielectric functions are composed of interband-transition-term (bound electron term) plus intraband-transition-term (free electron term: Drude term)

    2. Although you can obtain inerband-transition term from VASP, you should calculate intraband-transition-term by yourself especially for metallic system.

    3. It means that you need plasma frequency from the band structure.

    4. There are two possible method to calculate the plasma frequency using vasp code: see #14.

12. Dielectric constant for dielectrics with clamped-ion (update: 2016-01-28)

    1. Using the LEPSILON-tag, one can calculate the ion-clamped dielectric constant as written in VASP manual [LINK].

    2. Optimizing structure

    3. Calculating wavefunction

    4. Turning on LEPSILON-tag (LEPSILON = TRUE) and commenting out NPAR-tag or no NPAR tag in the INCAR-file.

13. Dielectric constant for dielectrics with relaxed-ion (update: 2016-01-28)

    1. Using the LEPSILON-tag & IBRION-tag, one can calculate the ion-clamped dielectric constant as written in VASP manual [LINK].

    2. Optimizing structure. Note that ionic forces should be very small because of the force constant calculations.

    3. Calculating wavefunction

    4. LEPSILON = TRUE / no NPAR tag / IBRION = 5,6,7,8 (one of them)

14. Prediction of plasma frequency (update: 2013-10-27)

    1. The plasma frequency I mention here is not the 'Screened plasma frequency'.

       1. You should distinguish between just plasma frequency and screened plasma frequency.

       2. The screened plasma frequency is the frequency where the real part of the dielectric function (inter-band term + intra-band term) becomes zero.

       3. And the screened plasma frequency is smaller than the plasma frequency.

    2. Using BoltzTraP code, you can obtain electrical conductivity over electron scattering time (SIGMA/TAU) within constant relaxation time approximation.

       1. Then, SIGMA/TAU can be transformed to plasma frequency.

       2. You shoud be careful when calculating plasma frequency.

       3. As the plasma frequency is obtaiend from the intra-band transition, I mean the group velocity of the electron,

       4. the k-point sampling is very important. You need sufficiently many k-points to obtain correct band dispersion.

       5. It is recommended to test the energy convergence with increasing number of k-points.

    3. 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.

15. Electrical conductivity using VASP code (update: 2013-11-13, 2015-04-17)

    1. Within semiclassical theory, one may obtain the electrical conductivity from the band structure by calculating the intraband structure.

    2. However, when you use the semiclassical theory, you should get the relaxation time from outside computations or you should calculate the relaxation time.

    3. If you already know the electron scattering time or relaxation time, you can get the electrical conductivity using VASP.

    4. There are two ways to get the electrical conductivity: see#14.

       1. One is using "BoltzTraP" code. And the other thing is to read the OUTCAR file from VASP code.

       2. 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.

       3. Unit of electrical conductivity per relaxation time is 10^22 [S/m/sec].

       4. 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

       5. 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. // 2015-04-17

       6. And the Energy dependent conductivity is also written in the "vasprun.xml" file. // 2015-04-17

16. MAGMOM tag in INCAR

    1. You may reduce the sentence length of MAGMOM tag as shown below,

       1. MAGMOM = 1 1 1 1 1 ==> MAGMOM = 5*1

       2. MAGMOM = 1 1 -3 -3 4 ==> MAGMOM = 2*1 2*-3 4

    2. It really saves your time to make INCAR file for spin polarized calculation of large supercell.

17. Many body perturbation theory (MBPT) : GW and BSE in VASP (I tested MBPT calculations with vasp version of 5.2.xx)

    1. MBPT calculations are very time consuming and very memory consuming.

    2. BSE tips in other site : http://cms.mpi.univie.ac.at/vasp-forum/forum_viewtopic.php?4.12605 (site broken)

    3. Calculation flow:

       1. DFT(Exc=LDA or GGA) calculations (or DFT+U, Hybrid-DFT) 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.

    4. Related tags:

       1. LOPTICS = T for derivative of wavefunction (see manual)

       2. NBANDS number of total bands including occupied and unoccupied, convergence test is needed (see manual)

       3. ENCUT convergence test is needed (see manual)

       4. CSHIFT smearing value for optical response results (freq. dependent dielectric function) (see manual)

       5. 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)

       6. NOMEGA (see manual)

       7. NBANDSO number of occupied bands you considered for BSE calculations

       8. NBANDSV number of unoccipied (virtual) bands you considered for BSE calculations

18. Charge analysis (update: 2015-04-20)

    1. Bader Charge

    2. You may analyze the atomic charge density distribution by using Bader Analysis. [http://theory.cm.utexas.edu/henkelman/code/bader/]

    3. The concept of Bader Charge is based on the space separation based on the minimum charge density positions.

    4. CHGCAR or CHG file

    5. You can directly analyze the atomic charge density using CHGCAR file.

    6. Or as mentioned in above (#6), you can get the atomic charge using DOSCAR or OUTCAR.

19. Unit conversion for VASP phonon frequency from density functional perturbation theory (DFPT) results, Unit of phonopy (update: 2015-06-24)

    1. [VASP]

       1. When one run DFPT calculations or BSE calculations, one obtain the vibration frequency in OUTCAR or vasprun.xml files.VASP.

       2. 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

    2. [phonopy] the default unit for frequency is THz. Therefore you can obtain phonon energy by producting planckconstant (h)

       1. 1 THz = 4.135669 meV

       2. example of PbTe phonon (VASP + phonopy, cubic lattice, lattice parameter of a=6.567843.

20. First principles calculations of Lattice thermal conductivity

    1. Thermal conductivity (kappa_tot) is a sum of electrical term and lattice term.

    2. 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)

    3. Lattice therm (kappa_latt) = (kappa_phonon) can be obtained by performing phono3py calculations, which calculates the transition rate of phonon-phonon scattering. (phono3py LINK)

21. Using phono3py code with VASP (2015-12-21)

    1. prepare the optimized supercell with the VASP parameter to be used in force calculations.

       1. The force should be less than 1.0e-08 eV/angstrom. Of course, you should optimized the structure considering the k-point mesh you will use in (1).

       2. WAVECAR file is needed to reduce computational resources for step (2).

    2. generate POSCAR files with pair-distance configurations. There might be one to ten thousands POSCAR files.

       1. phono3py -d --dim="4 3 3" -c POSCAR-unitcell --cutoff_pair_distance="4.4"

       2. supercell size, unitcell file, cutoff pair-distance in angstrom

    3. calculate force for pair-distance configurations using VASP.

       1. Make directories for configurations of POSCAR-00365.

       2. Copy INCAR, KPOINTS, POTCAR POSCAR, and WAVECAR.

       3. Using WAVECAR preapared from (0) will be useful to reduce the computational cost.

       4. In my case, I used the gamma point only but with the large supercell containing 200-400 atoms.

       5. Energy cutoff is also known to be sensitive to the lattice thermal conductivity tensor.

       6. For example, when the monoclinic unitcell volume is optimized within 400 eV and internal coordinate is optimized within 300 eV.

       7. then there could be non-zero off-digonal term in lattice thermal conductivity tensor.

    4. collect vasprun.xml file

       1. phono3py --cf3 disp-{00001..00999}/vasprun.xml (or) phono3py --cf3 disp-*/vasprun.xml

       2. disp_fc3.yaml, POSCAR-unitcell files should be in same directory.

    5. create fc2.hdf and fc3.hdf

       1. phono3py --dim="4 3 3" -c POSCAR-unitcell

       2. "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"

       3. The file size of fc3.hdf is very large. It is about 1 to 10 GB.

    6. calculate thermal conductivity

       1. split thermal 3phonon process configurations

          1. phono3py --fc3 --fc2 --dim="4 3 3" -v --mesh="11 11 11" -c POSCAR-unitcell --br --thm ---wgp

          2. q-mesh grid, write grid point for split calculation

          3. The computational time for 3phonon-transition rates are very time consuming when q-mesh grid is fine.

          4. Thermal conductivity, especially for low temperature region, q-mesh grid affects the lattice thermal conductivity.

          5. Therefore it is recommended to split the calculation set and calculate the transition rate separately using many cpu cores.

          6. After run the commands in (5.1), you will get "grid_address-m111111" and "ir_grid_points.yaml"

       2. calculate 3phonon transition rate for each q grid point

          1. phono3py --fc3 --fc2 --dim="4 3 3" -v --mesh="11 11 11" -c POSCAR-unitcell --br --thm --write_gamma --gp="3"

          2. Check q grid point number (gp) from ir_grid_points.yaml and write correct gp.

          3. The calculation time is dependent on q-mesh size.

          4. Each gp-point calculation time is varying from 10 minutes to 10 days depending supercell size and q-mesh size.

          5. Note that there are many gp if q-mesh size is fine.

       3. Sum up thermal conductivity

          1. phono3py --fc3 --fc2 --dim="4 3 3" -v --mesh="11 11 11" -c POSCAR-unitcell --br --thm --read_gamma > z_TC_def.txt

          2. Now you will get the lattice thermal conductivitiy file, "z_TC_def.txt"

          3. Type "tail -n 150 z_TC_def.txt" and find the lattice conductivity tensor as a function of temperature.

          4. (*) Useful tag: --isotope, --mv, --bmfp for isotope effect, mass variance effect, boundary scattering effect.

          5. (**) Auxilary tools

          6. kaccum and gaccum commands for accumulated lattice thermal conductivity and 3ph transition rate, respectively.

22. Band structure using VASP+Wannier90 (2016-04-05)

    1. 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

    2. Simple method

       1. Preconverge wavefunction

       2. Turn on LWANNIER90=T tag and run wan90 compiled VASP.

          1. Then Maximally Localized Wannier Functions (MLWFs) are generated with wannier90.win, wannier90.mmn, wannier90.eig, wannier90.out files.

          2. 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. (2017-03-24)

       3. Re-Run wannier90.

          1. wannier90.x wannier90

       4. Interpolate wannier functionals: first modify the wannier90.win file. For details, see the tutorial file.

          1. wannier90.x wannier90

    3. GW band structure from wannier90

       1. 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.

       2. 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).

       3. And you may know that NBANDS is dependent on NPAR tag. First check NBANDS without NPAR tag.

       4. Also you don't need too many wannier functions.

       5. 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

       6. 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).

       7. In shortly, you need to determine "NBANDS" for VASP and "NUM_WAN" for wannier90.

       8. Here I will use NBANDS=96 (I have 12 core nodes, so NBANDS should be divived by 12)

       9. Here I will use NUM_WAN = 8 (the number of wannier functions are related to the symmetries)

       10. Procedure for wannier

           1. (1-1) DFT calculations with default NBANDS

           2. (1-2) Optical calculations (LOPTICS=TRUE, NBANDS=96 )

           3. (1-3) GW calculations with wannier90 (ALGO=GW0, LWANNIER90=TRUE, no NPAR related tag)

              1. If you do not set wannier90.win file, then wannier90.win will be generated.

              2. But it is too heavy to play with it since there are too many BANDs.

              3. Make wannier90.win file (or modify wannier90.win file) as follows.

              4. (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)

              5. num_wann = 8

              6. num_bands = 8

              7. exclude_bands : 9-96

           4. (1-4) Run wannier90

              1. After getting wannier90.win, wannier90.eig, wannier90.mmn, rerun the wannier90.x

              2. cmd> wannier90.x wannier90

              3. Now you wannierize and obtain MLWFs.

           5. (1-5) Intepolate

              1. First modify or generate or make suitable wannier90.win or NAME.win.

              2. And run by typing "wannier90.x wannier90" or "wannier90.x NAME".

              3. Now you can interpolate for DOS, BANDS, or Fermi-Surface.

           6. IMPORTANT TIP for wannier90 with VASP (2016-11-02): 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.

23. Phonon thermal conductivity of low dimensional structure (2016/04/25)

    1. 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.

    2. 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.

    3. For example, I tested 1D carbon (C) atomic chain along z direction and the C-C distance is set to 1.28 angstrom.

    4. I considered 1x1x10 supercells containing 10 C atoms with various a=b values, 10, 14, 20, and 25 angstroms.

    5. 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.

    6. It implies that the cross section area should be checked carefully because the phono3py code assumes that one's structure is bulk.

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

    8. 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.