Using VASP for density functional theory (DFT)

Handy Terminal Commands

The commands you'll use every day on the server. Keys to press look like Esc; commands to type look like ls. Anything in [brackets] is a placeholder.

Connecting and moving files

Run scp in a terminal on your own computer, not while logged into the server.

Navigating and managing files

cp and mv overwrite files without asking, and rm is permanent (there's no trash). Names are case-sensitive: VASP reads INCAR, not incar.

Reading files

Running jobs and monitoring the server

Editing files with vim

vim opens in normal mode, where keys are commands. Press i for insert mode to type or paste, and Esc to go back. If keys start doing strange things, press Esc. Commands starting with : are typed in normal mode and appear at the bottom of the screen; press Enter to run them.

Moving around and editing (normal mode)

"Found a swap file" when opening a file means it's open somewhere else, or vim closed badly last time. If nobody else is editing it, choose E (Edit anyway) or D (Delete it).

General terminal tips

Before You Start: Creating the Necessary Files

Setting up your calculation folder

Every calculation gets its own folder holding five files: POSCAR, POTCAR, INCAR, KPOINTS and job. VASP reads them from the folder you submit the job from.

1. Log in and move into your working directory, e.g. cd DFT. Right after logging in you are in your home folder, one level too high. Run pwd to check where you are.
2. Make a folder for the compound and move into it, e.g. mkdir DyRuGa, then cd DyRuGa.
3. Create each file inside that folder with vim [filename]: press i, paste, then Esc + :wq. Don't paste file contents straight into the terminal prompt; the shell will try to run each line as a command.
4. I name the first folder as [composition-structure] and will have the specific composition named as composition or structure type, depending on what I’m doing, such as CeRuSi2-CeNiSi2. However, everyone works their own way, just make sure to keep detailed notes.

POSCAR file

1. Open CIF files with VESTA, which you can download here.
2. Export the structure in VASP format with fractional coordinates and name the file POSCAR.
   • On a Mac, VESTA saves it as POSCAR.vasp. Remove the extension with mv POSCAR.vasp POSCAR
3. Move POSCAR into your calculation folder on the server with scp (see Handy Terminal Commands).
4. Note the order of the elements on line 6 of POSCAR (e.g. Dy Ru Ga). The POTCAR must follow the same order.

POTCAR file

1. List the potentials available for an element with potshow, e.g. potshow Dy. Most elements have several versions, marked by a suffix (e.g. Dy and Dy_3; _pv and _sv also treat semicore states as valence). The VASP list of available PAW potentials explains which electrons each one treats as valence.
2. Use the potentials chosen for your project. For example, for DyRuGa we use Dy_3, Ru and Ga.
   • Dy_3 vs Dy: the _3 lanthanide potentials treat the atom as trivalent and freeze the 4f electrons in the core. This is cheaper and converges more easily, but the f electrons aren't described explicitly, so the results are more qualitative.
3. Build the POTCAR with potcreate, listing each element's potential once, in the same order as POSCAR. Only the potential names go in, not the number of atoms:
   • POSCAR elements Dy Ru Ga → potcreate Dy_3 Ru Ga
4. You can also generate POTCAR files using VASPKIT:
   • Type vaspkit in the terminal
   • Type 01 for VASP Input-Files Generator
   • Type 103 to generate POTCAR file with default settings
   • Note: Your PBE directory must be set up correctly when installing VASP
5. Check the result with grep TITEL POTCAR: the potentials should be listed in the same order as the elements in POSCAR.

INCAR file

The INCAR tells VASP which calculation to run. The template below is split into numbered sections. VASP ignores anything after a #, so you pick a calculation by uncommenting (deleting the # in front of) the tags in the sections you need, and leaving everything else commented.

1. [1] General Parameters is always active.
2. Activate one task section, [2]–[9], and comment out the others. [2] Structure Optimization is active in the template, so for any other task you must comment it out. With [2] off, VASP does a single static run (NSW defaults to 0).
3. The add-ons [10]–[13] (spin polarization, hybrid functionals, spin–orbit coupling, DFT+U) go on top of a task section when your system needs them.
4. Never leave the same tag active twice (e.g. two ISMEAR lines). Comment out the copy you aren't using.
5. Each calculation section on this page starts with a blue INCAR: note and the ready-to-use INCAR for that calculation ([1] plus the relevant sections). The table below summarizes them.

Copy the template below, or download it here. As written, it runs a geometry optimization ([1] + [2]).

### INCAR template for VASP calculations (adapted from Martin, 2015-12-22)
###
### HOW TO USE THIS TEMPLATE
###  - [1] General Parameters is ALWAYS active.
###  - Activate ONE task section [2]-[9] by deleting the # in front of its tags,
###    and comment out (#) every other task section, including [2].
###    With [2] commented out, NSW defaults to 0 (a static run).
###  - Add-ons [10]-[13] can be combined with a task section when needed.
###  - Never leave the same tag active twice (e.g. two ISMEAR lines).
###    Comment out the copy you are not using.
###  - As written, this file runs a geometry optimization ([1] + [2]).

## [1] General Parameters / Electronic Relaxation -- ALWAYS ACTIVE
PREC   = Accurate  # precision
ENCUT  = 500       # plane-wave cutoff energy (eV)
EDIFF  = 1E-8      # convergence criterion for the electronic loop
NELM   = 150       # maximum number of electronic iterations
NELMIN = 4         # minimum number of electronic iterations
KPAR   = 2         # k-points run in parallel (comment out for Bader)
NCORE  = 4         # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [2] Structure Optimization -- ACTIVE BY DEFAULT (comment out for any other task)
EDIFFG = 1E-6      # convergence criterion for the ionic loop (pos = energy; neg = forces)
NSW    = 1000      # maximum number of ionic steps
IBRION = 1         # ionic update: 1 = quasi-Newton, 2 = conjugate gradient (for bad guesses)
ISIF   = 3         # freedom: 2 = ions; 4 = ions + shape; 3 = ions + shape + volume
MAXMIX = 80        # max steps stored in Broyden mixer (negative: resets after each ionic step)
ISMEAR = 1         # smearing: -5 = insulators/semiconductors/DOS; 0 = large systems (SIGMA = 0.05); 1 = metals (SIGMA = 0.2)

## [3] DOS and Band Structure
#NEDOS  = 1001     # number of DOS grid points
#LORBIT = 10       # write DOSCAR and PROCAR (11 for lm-decomposed orbitals)
#ISMEAR = -5       # tetrahedron method for DOS
# if you add ICHARG = 11, a previous CHGCAR must be in the folder

## [4] Bader Charges and Electron Localization Function
#LAECHG = .TRUE.   # write AECCAR0 + AECCAR2 (core + valence charge) for Bader
#LELF   = .TRUE.   # write ELFCAR
#NPAR   = 1

## [5] COHP Analysis (via LOBSTER)
#LWAVE  = .TRUE.   # LOBSTER needs the WAVECAR
#NSW    = 0        # static run
#ISYM   = -1       # symmetry off
#NBANDS = 600      # should equal the total number of orbitals (see Troubleshooting)
#LORBIT = 12
#ISTART = 0

## [6] Phonon Dispersion via finite displacements
#EDIFF   = 1E-8    # replaces EDIFF in [1] -- don't leave both active
#ISMEAR  = 0
#SIGMA   = 0.01
#ADDGRID = .TRUE.

## [7] Stress Tensor
#ENCUT  = 600      # replaces ENCUT in [1] -- comment that one out
#IBRION = 6
#ISIF   = 6
#NSW    = 1
#NFREE  = 2
## additionally for DFPT (replaces IBRION = 6 above):
#IBRION = 8
#NWRITE = 3

## [8] Supercell (band offset) calculations
#LVHAR = .TRUE.
#AMIN  = 0.1

## [9] Born Effective Charge
#LEPSILON = .TRUE. # dielectric matrix
#LOPTICS  = .TRUE. # DFPT including empty bands for hybrid functionals

## ---------- ADD-ONS (combine with a task section above) ----------

## [10] Spin Polarization / Magnetism
#ISPIN  = 2        # spin-polarized calculation
#MAGMOM = 4*0 4*1  # starting moment per atom, in POSCAR order (here 4 atoms at 0, then 4 at 1)

## [11] Hybrid Calculations -- ONLY add to [2] if a normal optimization won't converge (very slow)
#LHFCALC  = .TRUE. # use HF/DFT hybrid functionals
#AEXX     = 0.25   # fraction of exact exchange
#AGGAX    = 0.75   # fraction of GGA exchange
#ALGO     = Damped # electronic minimization algorithm (damped velocity friction)
#TIME     = 0.05   # time step for ALGO = Damped (lower it if convergence fails)
#PRECFOCK = FAST   # FFT grid for the exact exchange
#HFSCREEN = 0.2    # range-separation parameter
#NKREDX   = 2      # grid reduction factors: each KPOINTS number must be divisible by these
#NKREDY   = 2
#NKREDZ   = 2

## [12] Spin-Orbit Coupling -- run with vasp_ncl, not vasp_std
#LSORBIT       = .TRUE.
#LNONCOLLINEAR = .TRUE.

## [13] GGA+U, Coulomb interaction energy (see Hubbard U Correction)
#LDAU      = .TRUE.
#LDAUTYPE  = 1     # 1 = Liechtenstein et al.; 2 = Dudarev (only U - J is meaningful)
#LDAUPRINT = 0
#LDAUL     = -1 2 -1 -1  # orbital to correct, one per species in POSCAR order: -1 = none, 2 = d, 3 = f
#LDAUU     = 0 3.71 0 0  # U (eV), one per species
#LDAUJ     = 0 0 0 0     # J (eV), one per species
#LASPH     = .TRUE.

KPOINTS file

KPOINTS sets how finely VASP samples reciprocal space (the Brillouin zone). A denser mesh (bigger numbers) gives more accurate results but takes longer. Longer lattice vectors need fewer k-points along that direction.

1. Generate the file with VASPKIT: type vaspkit, then 01 (VASP Input-Files Generator), then 102 (KPOINTS). When asked for the k-point resolution, use the Fine level, 0.02. See the geometry optimization video for a walkthrough.
2. Inspect the generated file with cat KPOINTS (or vim KPOINTS). It looks like this:

   K-Points
   0
   Gamma
   8 2 8
   0 0 0

   • Line 1 (K-Points): a comment; any text works.
   • Line 2 (0): 0 tells VASP to generate the mesh automatically.
   • Line 3 (Gamma): the mesh is centered on the Γ point, the center of the Brillouin zone (Monkhorst-Pack is the alternative).
   • Line 4 (8 2 8): the number of k-points along each reciprocal lattice vector: 8 along the first, 2 along the second and 8 along the third.
   • Line 5 (0 0 0): an optional shift of the mesh; here none.
3. Make sure all three mesh numbers on line 4 are even. We run with KPAR = 2, which splits the k-points into two groups, so an even mesh divides evenly between them. If VASPKIT gives an odd number, round it up in vim (e.g. 7 3 7 → 8 4 8). Even numbers are also divisible by the NKRED factors used in hybrid calculations ([11]; see Troubleshooting).

job file (SLURM)

1. Copy the job template into your calculation folder, or copy one from a previous calculation: cp ../[otherfolder]/job .
2. Edit it with vim job:
   • #SBATCH -J: give each run its own job name so it's easy to find in sq and cancel.
   • #SBATCH -n: the number of cores. Use the number you've been given (e.g. 16), not the template's 32.
3. #!/bin/sh must be the very first line of the file. Copy-pasting often adds a blank or stray line above it, and sbatch then refuses the file with an error like "This does not look like a batch script. The first line must start with #!". Open the file with vim job, put the cursor on the extra line, press dd, then :wq.
4. VASP's output goes to a file called out (from the last line, srun vasp_std > out). After submitting with sbatch job, follow it with tail -f out.

Final check before submitting

• ls in your folder shows INCAR  job  KPOINTS  POSCAR  POTCAR.
• The POTCAR order matches POSCAR line 6 (grep TITEL POTCAR).
• Only the INCAR sections for your calculation are active. [11] Hybrid Calculations stays commented out unless you need it.
• All KPOINTS mesh numbers are even.
• If you're new, have someone look over your files before you run.

Installing VASPKIT

1. Visit the VASPKIT installation page for step-by-step instructions for Linux, macOS, and Windows.
2. Follow the guide to download, compile, and set up VASPKIT on your system.
3. For a visual walkthrough, see the tutorial video above.
4. Once installed, you can run VASPKIT by typing vaspkit in your terminal.

Geometry Optimization

INCAR: keep [1] General Parameters and [2] Structure Optimization active (the template as written). Every other section stays commented out; in particular, make sure [11] Hybrid Calculations is commented out.

# INCAR for geometry optimization: [1] + [2]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [2] Structure Optimization
EDIFFG   = 1E-6        # convergence criterion for the ionic loop (pos = energy; neg = forces)
NSW      = 1000        # maximum number of ionic steps
IBRION   = 1           # ionic update: 1 = quasi-Newton, 2 = conjugate gradient (for bad guesses)
ISIF     = 3           # freedom: 2 = ions; 4 = ions + shape; 3 = ions + shape + volume
MAXMIX   = 80          # max steps stored in Broyden mixer (negative: resets after each ionic step)
ISMEAR   = 1           # smearing: -5 = insulators/semiconductors/DOS; 0 = large systems (SIGMA = 0.05); 1 = metals (SIGMA = 0.2)

1. Note: mpirun is a command used to launch parallel programs written using the Message Passing Interface (MPI) standard.
2. The Oliynyk lab uses SLURM for job scheduling, so you need a job file. See the template here and the job file notes above.
3. Once all the above files are created, you’re ready to go with optimizing. You MUST do this before anything else. It is pertinent to start any calculations with an optimized and most thermodynamically stable geometry.
4. To run with SLURM, use:
   • sbatch job
   • Check status: sq
   • Cancel a job: scancel [job no.] (job number shown by sq)
   • Tip: Change the job name in your job file for each run so cancelling is easy.
5. type mpirun -np 16 vasp_std>1& (Oliynyk lab - use sbatch job after setting job file) in the command line while in the directory of the folder with your files. "1" indicates the first run, "2" indicates the second, and "3" indicates the third. I do this to keep track of everything, but you can just repeat the same command after copying CONTCAR to POSCAR.
6. Can view status of calculation/warnings with the following command:
   • Vim + 1
   • tail 1 (only displays last couple of lines)
   • With sbatch job, the output file is out: tail -f out
7. Can view energy with the following command:
   • grep TOTEN OUTCAR
8. If ‘reached required accuracy - stopping structural energy minimisation’ is displayed at the bottom of file 1, this indicates it has reached your threshold of accuracy, set by EDIFFG in [2] of the INCAR file. Copy the optimized structure to the structure input file with:
   • cp CONTCAR POSCAR
9. Continue with the second run:
   • mpirun -np 4 vasp_std>2& (Oliynyk lab - use sbatch job after setting job file)
10. If ‘reached required accuracy - stopping structural energy minimisation’ is displayed at the bottom of file 2, copy the optimized structure to the structure input file with:
    • cp CONTCAR POSCAR
11. Continue with the third run:
    • mpirun -np 4 vasp_std>3& (Oliynyk lab - use sbatch job after setting job file)
12. If ‘reached required accuracy - stopping structural energy minimisation’ is displayed at the bottom of file 3, copy the optimized structure to the structure input file with:
    • cp CONTCAR POSCAR

For hybrid calculations

INCAR: [1] General Parameters + [2] Structure Optimization + [11] Hybrid Calculations.

# INCAR for hybrid optimization: [1] + [2] + [11]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [2] Structure Optimization
EDIFFG   = 1E-6        # convergence criterion for the ionic loop (pos = energy; neg = forces)
NSW      = 1000        # maximum number of ionic steps
IBRION   = 1           # ionic update: 1 = quasi-Newton, 2 = conjugate gradient (for bad guesses)
ISIF     = 3           # freedom: 2 = ions; 4 = ions + shape; 3 = ions + shape + volume
MAXMIX   = 80          # max steps stored in Broyden mixer (negative: resets after each ionic step)
ISMEAR   = 1           # smearing: -5 = insulators/semiconductors/DOS; 0 = large systems (SIGMA = 0.05); 1 = metals (SIGMA = 0.2)

## [11] Hybrid Calculations
LHFCALC  = .TRUE.      # use HF/DFT hybrid functionals
AEXX     = 0.25        # fraction of exact exchange
AGGAX    = 0.75        # fraction of GGA exchange
ALGO     = Damped      # electronic minimization algorithm (damped velocity friction)
TIME     = 0.05        # time step for ALGO = Damped (lower it if convergence fails)
PRECFOCK = FAST        # FFT grid for the exact exchange
HFSCREEN = 0.2         # range-separation parameter
NKREDX   = 2           # grid reduction factors: each KPOINTS number must be divisible by these
NKREDY   = 2
NKREDZ   = 2

1. If a certain potential cannot converge no matter what, uncomment [11] Hybrid Calculations on top of [1] + [2]. The optimization will take much longer, up to weeks in some cases.
2. Each KPOINTS mesh number must be divisible by NKREDX, NKREDY and NKREDZ (2), i.e. even.
3. In case there’s a problem with convergence, lower TIME in [11], e.g. halve it from 0.05 to 0.025.
4. Setting an easier convergence criterion may help with stubborn structures: EDIFF = 1E-6 in [1] (electronic loop) and EDIFFG = 1E-5 in [2] (ionic loop).

Data Accrual

1. For free energy:
   • While in the structure directory, type grep TOTEN OUTCAR to see the final energy of the optimized structure, which is the last iterated value at the bottom.
   • Looks like this:
   
2. For optimized parameters:
   • While in the structure directory, type vim POSCAR
   • This is why we type cp CONTCAR POSCAR at the end of our third optimization iteration.
   • One can also just type vim CONTCAR, but we want to export POSCAR later for further calculations – so it’s best to do the first command in 1.
   • Looks like this:
   
   • After opening the file in the vim editor, you will see the optimized parameters as displayed in the image above.
   • Where the first line with values indicates parameter A, second indicates B, and third indicates C (in Angstroms).
   • In this case, we have a matrix representation of a non-orthogonal unit cell. It is monoclinic. You will notice that there are more than one value in each line.
   
   • To deal with this and convert it to appropriate parameters as well as angles, import the optimized POSCAR to VESTA and report the A, B and C cell parameters as well as the β Angle (°) in the summary section.
     

Density of States

INCAR: [1] General Parameters + [3] DOS and Band Structure. Comment out [2] Structure Optimization. Keep LWAVE = .FALSE. commented, as we need the WAVECAR file.

# INCAR for density of states: [1] + [3]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [3] DOS and Band Structure
NEDOS    = 1001        # number of DOS grid points
LORBIT   = 10          # write DOSCAR and PROCAR (11 for lm-decomposed orbitals)
ISMEAR   = -5          # tetrahedron method for DOS
# if you add ICHARG = 11, a previous CHGCAR must be in the folder

1. Non-spin polarized DOSCAR
2. Spin-Polarized and Non-Spin-Polarized :
   • ISPIN=1 (non-spin-polarized) and ISPIN=2 (spin-polarized) calculations.
3. Atomic orbital contribution settings:
• LORBIT=11-14: Individual orbitals are lm-resolved/decomposed (e.g., px, py, pz, etc.).
• LORBIT=0, 1, 2, 5, 10: Orbitals are grouped (e.g., total p-orbital contributions).
4. Run with mpirun -np 16 vasp_std>Dos1& (Oliynyk lab - use sbatch job after setting job file)
5. Copy DOSCAR files to local desktop for analysis when converged
6. Plotting DOSCAR:
   • You can plot the exported POSCAR and DOSCAR files using DOSCAR-GUI.
   • More information and settings can be found on GitHub.

Magnetic Moment Calculation

1. Spin-polarized DOSCAR (only if magnetism is expected here, otherwise you can ignore this)
   • INCAR: the DOS setup ([1] + [3]) plus [10] Spin Polarization:

   # INCAR for spin-polarized DOS: [1] + [3] + [10]
   # (all other template sections commented out or deleted)
   
   ## [1] General Parameters
   PREC     = Accurate    # precision
   ENCUT    = 500         # plane-wave cutoff energy (eV)
   EDIFF    = 1E-8        # convergence criterion for the electronic loop
   NELM     = 150         # maximum number of electronic iterations
   NELMIN   = 4           # minimum number of electronic iterations
   KPAR     = 2           # k-points run in parallel (comment out for Bader)
   NCORE    = 4           # cores working on each orbital
   #LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP
   
   ## [3] DOS and Band Structure
   NEDOS    = 1001        # number of DOS grid points
   LORBIT   = 10          # write DOSCAR and PROCAR (11 for lm-decomposed orbitals)
   ISMEAR   = -5          # tetrahedron method for DOS
   # if you add ICHARG = 11, a previous CHGCAR must be in the folder
   
   ## [10] Spin Polarization / Magnetism
   ISPIN    = 2           # spin-polarized calculation
   MAGMOM   = 4*0 4*1     # starting moment per atom, in POSCAR order (here 4 atoms at 0, then 4 at 1)

   • ISPIN = 2 turns on spin polarization
   • MAGMOM gives a starting moment for each atom, in POSCAR order: 4*0 4*1 is the same as 0 0 0 0 1 1 1 1
   • Whole numbers are fine as a starting guess
   • This is your prediction, the OUTCAR file will give you the software's guess
   • Run with mpirun -np 8 vasp_std>Dos2& (Oliynyk lab - use sbatch job after setting job file)
   • Rerun until the OUTCAR charges are the same as your input MAGMOM prediction, as it gets better with better predictions
   • As in they’re all positive numbers or negative numbers with respect to your entry. The output does not have to match besides it being positive or negative, usually.
   • Retrieve DOSCAR file when converged

Electron Localization Function

INCAR: [1] General Parameters + [4] Bader Charges and ELF. Comment out [2] Structure Optimization and KPAR in [1].

# INCAR for ELF and Bader: [1] + [4]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
#KPAR     = 2          # commented out for Bader
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [4] Bader Charges and Electron Localization Function
LAECHG   = .TRUE.      # write AECCAR0 + AECCAR2 (core + valence charge) for Bader
LELF     = .TRUE.      # write ELFCAR
NPAR     = 1

ELF and Bader use the same run. [4] writes the ELFCAR (from LELF) and the AECCAR0/AECCAR2 files Bader needs (from LAECHG) at the same time, so there's no need to run VASP twice. Run it once, then do whichever post-processing you need: open ELFCAR for the ELF (below), or follow the Bader Charges steps to get ACF.dat.

1. Run with mpirun -np 16 vasp_std>Elf& (Oliynyk lab - use sbatch job after setting job file)
2. No extra commands are needed for the ELF: import the ELFCAR file into your local computer and drag it into VESTA to view

Bader Charges

INCAR: same as ELF: [1] General Parameters + [4] Bader Charges and ELF, with [2] Structure Optimization and KPAR commented out (see the INCAR in Electron Localization Function). If you already did the ELF run, skip straight to chgsum.pl below in that folder.

1. Run with mpirun -np 8 vasp_std>Bader& (Oliynyk lab - use sbatch job after setting job file)
2. Unlike the ELF, Bader needs two post-processing commands to produce ACF.dat. When converged, run chgsum.pl AECCAR0 AECCAR2
   • This combines core and valence electrons
3. type bader CHGCAR -ref CHGCAR_sum, which writes ACF.dat
   • Looks like this:
   
4. Import ACF.dat to local machine or just view it in the terminal with vim command
5. Run grep ZVAL POTCAR to get the number of valence electrons (ZVAL) of each element
   • Looks like this:
   
   • ZVAL values are listed in the POTCAR order, which is the same as the element order in POSCAR
6. Match each row of ACF.dat to its atom. Rows are in the same order as the atoms in POSCAR (vim POSCAR if you forget). For example, if POSCAR lists Dy Ru Ga with counts 4 4 4, rows 1–4 are Dy, rows 5–8 are Ru and rows 9–12 are Ga.
7. For each atom, subtract the Bader CHARGE from that element's ZVAL: net charge = ZVAL − CHARGE. A positive result means the atom lost electrons (cation-like); a negative result means it gained them (anion-like).

Worked example

Row 1 of ACF.dat:

    #         X           Y           Z       CHARGE      MIN DIST   ATOMIC VOL
 --------------------------------------------------------------------------------
    1    1.8293      1.8293      1.9646      7.8550      1.3155     17.8112

Atom 1 is the first atom in POSCAR, Dy, and grep ZVAL POTCAR gives ZVAL = 9 for the Dy_3 potential. Read the CHARGE column (the fifth number), not Z:

So Dy has given away about 1.15 electrons. Repeat this for every row, using the ZVAL of that row's element. As a check, NUMBER OF ELECTRONS at the bottom of ACF.dat should equal the sum of ZVAL over all atoms.

Crystal Orbital Hamiltonian/Overlap Population

INCAR: [1] General Parameters + [5] COHP Analysis. Comment out [2] Structure Optimization. Set NBANDS for your system (see the note below and calculating NBANDS).

# INCAR for COHP/COOP: [1] + [5]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [5] COHP Analysis
LWAVE    = .TRUE.      # LOBSTER needs the WAVECAR
NSW      = 0           # static run
ISYM     = -1          # symmetry off
NBANDS   = 600         # should equal the total number of orbitals (see Troubleshooting)
LORBIT   = 12
ISTART   = 0

1. Copy the lobsterin file from this template below, or download it here. You can also add the lobsterin file after your calculation. lobsterin is not needed for the VASP calculation to run. Only for lobster to run.
! First, enter the energetic window in eV (relative to the Fermi level):
COHPstartEnergy -10
COHPendEnergy 5
!
! Then, specify which types of valence orbitals to use:
includeorbitals s p d
! You can also specify the basis functions per element manually, e.g.:
! basisfunctions Ga 4s 4p
! basisfunctions Sr 5s 4s 4p ! Sr_sv potential
!
! Now define the pairs for which COHP analysis etc. should be done.
! The atoms are numbered as per order in the PAW-code control file.
!cohpbetween atom 1 atom 10
!
! If you are interested in single orbital COHPs, you can get all the pairs ! like s-s, s-p_x, ..., p_z-p_z. Uncomment this line to switch it on:
! cohpbetween atom 1 atom 2 orbitalwise
!
! If you want to generate the COHPPairs automatically, use this to include ! all pairs in a given distance range (in Angstrom, not in atomic units):
cohpGenerator from 1 to 2.5 type Si type Fe
cohpGenerator from 1 to 3.3 type Gd type Si
cohpGenerator from 1 to 3.3 type Gd type Fe
cohpGenerator from 1 to 2.6 type Si type Si
! and in the latter case only between the specified elements
- of course adjust accordingly to your atomic interactions and optimized structure's bond distances.
- I like using VESTA to visualize this

2. type mpirun -np 20 vasp_std>COHP& (Oliynyk lab - use sbatch job after setting job file)
3. type lobster
4. Note: If you do not set NBANDS in [5] or set it too low, you may see an NBANDS error in the output file. Check the required value, then set NBANDS to that value plus 20 for reliable results.
5. Import COOPCAR.lobster, COHPCAR.lobster, to local machine as needed for plotting
6. To plot COHP/COOP, save the files (COOPCAR.lobster, COHPCAR.lobster, etc.) in a .zip archive and upload to COHP-COOP-GUI.
7. See GitHub for more info and settings.

Band Structure

INCAR: same as DOS: [1] General Parameters + [3] DOS and Band Structure. Comment out [2] Structure Optimization. Keep LWAVE = .FALSE. commented, as we need the WAVECAR file.

# INCAR for band structure: [1] + [3]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [3] DOS and Band Structure
NEDOS    = 1001        # number of DOS grid points
LORBIT   = 10          # write DOSCAR and PROCAR (11 for lm-decomposed orbitals)
ISMEAR   = -5          # tetrahedron method for DOS
# if you add ICHARG = 11, a previous CHGCAR must be in the folder

1. We need to use VASPKIT to create a new KPOINTS using the symmetry. Follow these steps:
• Type vaspkit into the terminal/powershell
• Type 03 to select K-Path for Band-Structure
• Type 303 for bulk structure
• Run cp KPOINTS KPOINTS_old to back up your old KPOINTS
• Run cp KPATH.in KPOINTS to use the new K-path
• This new KPATH.in will become your KPOINTS for the band structure only. VASPKIT automatically suggests a K-Path based on your structure's symmetry.

Effective Mass Calculation with VASPKIT

This section covers how to calculate effective masses for electrons and holes using VASPKIT, which is particularly useful for semiconductor materials and electronic transport properties.

1. Run Regular Relaxation

First, perform a standard geometry optimization as described in the Geometry Optimization section to obtain the optimized structure.

2. Run Band Structure Calculation (SCF)

Perform a band structure calculation following the Band Structure section instructions. This provides the electronic band structure needed for effective mass calculations.

3. Determine High-Symmetry Points of Interest

Identify the high-symmetry points in your material's band structure where you want to calculate effective masses, typically at band extrema (valence band maximum and conduction band minimum).

4. Use VASPKIT Option 911 for Band Edge Analysis

VASPKIT option 911 can automatically determine band edge positions, namely VBM (Valence Band Maximum) and CBM (Conduction Band Minimum):

+-------------------------- Warm Tips --------------------------+
ONLY Reliable for Band-Structure Calculations!
+---------------------------------------------------------------+
-->> (01) Reading Input Parameters From INCAR File...
-->> (02) Reading Energy-Levels From EIGENVAL File...
+-------------------------- Summary ----------------------------+
Band Character: Direct
Band Gap (eV): 1.7711
Eigenvalue of VBM (eV): -2.1236
Eigenvalue of CBM (eV): -0.3526
HOMO & LUMO Bands: 9 10
Location of VBM: 0.333333 0.333333 0.000000
Location of CBM: 0.333333 0.333333 0.000000
+---------------------------------------------------------------+
5. Prepare INPUT.in File for Effective Mass Calculation

Next, calculate the effective mass of electrons and holes along specific directions (e.g., K → Γ and K → M at point K). Create an INPUT.in file with the following content:

1 # Setting it to 1 will generate a KPOINTS file for calculating the effective mass; setting it to 2 will calculate the effective mass.
6 # 6 indicates that we take 6 discrete points near K for polynomial fitting.
0.015 # k-cutoff, cutoff radius, typical value 0.015 1/Å
3 # 3 indicates there are three effective mass tasks.

0.500000 0.500000 0.500000 0.000000 0.000000 0.500000 # A->Z Calculate effective mass from A toward Z.
0.500000 0.500000 0.500000 0.000000 0.500000 0.500000 # A->R Calculate effective mass from A toward R.
0.500000 0.500000 0.500000 0.500000 0.500000 0.000000 # A->M Calculate effective mass from A toward M.


Parameters explained:

• Line 1: Mode selector (1 = generate KPOINTS, 2 = calculate effective mass)
• Line 2: Number of discrete k-points for polynomial fitting (typically 6-8)
• Line 3: k-cutoff radius in 1/Å (typical value: 0.015)
• Line 4: Number of k-path directions to analyze
• Lines 5-7: Start and end points for each k-path in fractional coordinates
6. Generate KPOINTS and Run VASP

With the first line set to 1 in INPUT.in, run VASPKIT option 913 to generate the appropriate KPOINTS file, then run VASP with this new KPOINTS file. (It will overwrite KPOINTS in the VASP working directory, so save a backup of your old KPOINTS file.)

7. Calculate Effective Mass

After the VASP calculation is completed, modify the first line of INPUT.in from 1 to 2, then run VASPKIT option 913 again to calculate the effective masses:

Band Index: LUMO = 10 HOMO = 9
Effective-Mass (in m0) Electron (Prec.) Hole (Prec.)
K-PATH No.1: K->G 0.473 (0.3E-07) -0.569 (0.1E-07)
K-PATH No.2: K->M 0.524 (0.1E-06) -0.691 (0.1E-06)
8. Results Interpretation
• Effective masses are given in units of free electron mass (m₀)
• Electron effective mass refers to the conduction band minimum
• Hole effective mass refers to the valence band maximum
• Precision values in parentheses indicate the fitting accuracy
• Lower effective mass indicates higher carrier mobility
9. Important Notes
• Effective mass calculations require high-quality band structure data
• The k-cutoff radius should be small enough to capture local curvature but large enough for accurate fitting
• Results are only reliable near band extrema (VBM and CBM)
• For anisotropic materials, calculate along different crystallographic directions
• Always verify the band gap and band positions before calculating effective masses

Boltzmann Transport Properties with VASPKIT

Reference tutorial

This section covers how to calculate electronic transport properties using the Boltzmann transport equation with VASPKIT, including electrical conductivity, Seebeck coefficient, power factor, and carrier concentration.

1. Run Regular Relaxation

First, perform a standard geometry optimization as described in the Geometry Optimization section to obtain the optimized structure.

2. Run Band Structure Calculation (SCF)

Perform a band structure calculation following the Band Structure section instructions. This provides the electronic band structure needed for transport property calculations.

3. Generate KPOINTS for Transport Properties

Call the VASPKIT-681 command to generate a KPOINTS file specifically for the calculation of transport properties:

vaspkit
681

This will generate the appropriate k-point mesh needed for accurate transport property calculations.

4. Prepare INPUT.in File

Create an INPUT.in file (note that VASPKIT.in will be renamed to INPUT.in) with the following parameters for calculating transport properties:

1 # 1 for gradient method with constant relaxation time approximation
3D # 2D for slab, 3D for bulk
-2.0 2.0 # Minimum and maximum values of chemical potential with respect to Fermi energy
1000 # Number of intervals between the minimum and maximum chemical potential values
300.0 # Temperature (float type, in units of K)
1.0 # Relaxation time (float type, in units of s), Set 1.0 if you don't know the exact relaxation time
0.0 # Desired gap value in scissor correction (float type, in units of V), Zero means make no correction, available only non-spin-polarization calculation
1 # Finite difference method (integer type)

Parameters explained:

• Line 1: Calculation method (1 = gradient method with constant relaxation time)
• Line 2: System dimensionality (2D for slabs, 3D for bulk materials)
• Line 3: Chemical potential range in eV relative to Fermi energy
• Line 4: Number of energy intervals for integration
• Line 5: Temperature in Kelvin
• Line 6: Relaxation time in seconds (use 1.0 as default if unknown)
• Line 7: Scissor correction for band gap (0.0 = no correction)
• Line 8: Finite difference method selection
5. Calculate Transport Properties

Call the VASPKIT-682 command to calculate the transport properties:

vaspkit
682
-->> (01) Reading INPUT.in File...
-->> (02) Reading Structural Parameters from POSCAR File...
-->> (03) Reading Input Parameters From INCAR File...
+---------------------------------------------------------------+
| >>> The Fermi Energy will be set to zero eV <<< |
+---------------------------------------------------------------+
-->> (04) Reading Fermi-Energy from DOSCAR File...
-->> (05) Reading Energy-Levels From EIGENVAL File...
+-------------------------- Summary ----------------------------+
Chemical Potential Range (eV): -2.00 2.00
Number of Intervals: 1000
Temperature (K): 300.00
Relaxation Time (s): 0.10E+01
+---------------------------------------------------------------+
-->> (06) Calculating Boltzmann Transport Properties...
Percentage complete: 25.0%
Percentage complete: 50.0%
Percentage complete: 75.0%
Percentage complete: 100.0%
-->> (07) Written ELECTRONIC_CONDUCTIVITY.dat File!
-->> (08) Written SEEBECK_COEFFICIENT.dat File!
-->> (09) Written POWER_FACTOR.dat File!
-->> (10) Written ELECTRONIC_THERMAL_CONDUCTIVITY.dat File!
-->> (11) Written CARRIER_CONCENTRATION.dat File!
6. Output Files and Analysis

The calculation will generate several output files:

• ELECTRONIC_CONDUCTIVITY.dat - Electrical conductivity vs chemical potential
• SEEBECK_COEFFICIENT.dat - Seebeck coefficient vs chemical potential
• POWER_FACTOR.dat - Thermoelectric power factor vs chemical potential
• ELECTRONIC_THERMAL_CONDUCTIVITY.dat - Electronic thermal conductivity vs chemical potential
• CARRIER_CONCENTRATION.dat - Carrier concentration vs chemical potential
7. Data Visualization

You can plot the results using Python. Here's a sample code to plot carrier concentration:

import numpy as np
import matplotlib.pyplot as plt

# Read the data from CARRIER_CONCENTRATION.dat
data = np.loadtxt('CARRIER_CONCENTRATION.dat', skiprows=1) # Skip the header row

# Extract energy and carrier concentration
energy = data[:, 0] # Energy in eV
carrier_concentration = data[:, 1] # Carrier concentration in 10^20/cm³

# Create the plot
plt.figure(figsize=(10, 6))
plt.plot(energy, carrier_concentration, 'b-', linewidth=2, label='Carrier Concentration')
plt.xlabel('Energy (eV)', fontsize=12)
plt.ylabel('Carrier Concentration (10²⁰/cm³)', fontsize=12)
plt.title('Carrier Concentration vs Energy', fontsize=14, fontweight='bold')
plt.grid(True, alpha=0.3)
plt.legend()

# Add a horizontal line at y=0 for reference
plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)

# Improve layout
plt.tight_layout()

# Show the plot
plt.show()

# Optional: Save the plot
plt.savefig('carrier_concentration_plot.png', dpi=300, bbox_inches='tight')
print("Plot saved as 'carrier_concentration_plot.png'")

Example Output:


Example carrier concentration plot generated using the code above
8. Important Notes
• Transport properties are highly sensitive to the quality of the band structure calculation
• Use dense k-point meshes for accurate results
• The relaxation time parameter significantly affects absolute values but not trends
• For thermoelectric applications, focus on the power factor and ZT calculations
• Temperature-dependent calculations can reveal optimal operating conditions
• Always validate results against experimental data when available

Calculating Hubbard U Correction For DFT+U (LSDA+U)

1. File Setup for Selective U Application

Important: When applying Hubbard U corrections to only some atoms of the same element type (e.g., applying +U to only some U atoms in URhGe), you must modify both POSCAR and POTCAR files to distinguish between atoms that receive the correction and those that don't.

Example for URhGe system:

If your original POSCAR has:

URhGe
4.0
...
U  Rh  Ge
4  4   4
            

And you want to apply +U to only the first U atom, modify it to:

URhGe
4.0
...
U  U   Rh  Ge
1  3   4   4
            

This creates two distinct "U" species: the first one (1 atom) will receive the +U correction, and the second one (3 atoms) will be treated as regular U atoms.

Correspondingly, your POTCAR file must be created with:

potcreate U U Rh Ge

This ensures that VASP recognizes four distinct atomic species, allowing you to apply different parameters to each.

With this setup, your LDA+U parameters would be:

LDAU         = .TRUE.
LDAUTYPE     = 3
LDAUL        = 3  -1  -1  -1     # Apply to f-orbitals of 1st U only
LDAUU        = 3.42 0.00 0.00 0.00  # U parameter only for 1st U species
LDAUJ        = 0.51 0.00 0.00 0.00  # J parameter only for 1st U species
            

The first set of parameters applies to the first U species (receiving +U), while the remaining parameters apply to the second U species, Rh, and Ge respectively (all set to not receive corrections).

2. Ground State Calculation

Prepare the required input files for your system. The example below shows NiO with antiferromagnetic ordering, but the approach applies similarly to other systems like UGe₂.

POSCAR setup (following the file setup from step 1):

For NiO where you want to apply +U to all Ni atoms but distinguish different magnetic sites:

NiO AFM
4.03500000 
2.0000000000   1.0000000000   1.0000000000 
1.0000000000   2.0000000000   1.0000000000 
1.0000000000   1.0000000000   2.0000000000 
Ni  Ni  O
1   15   16 (adding one more Ni species to distinguish sites, so total 2 Ni species) 
Direct
...atomic positions...
            

POTCAR setup:

Create the POTCAR file with distinct entries for each species defined in POSCAR:

potcreate Ni Ni O

INCAR file for ground state calculation:

    SYSTEM       = NiO AFM
    PREC         = A
    EDIFF        = 1E-6
    ISMEAR       = 0
    SIGMA        = 0.2
    ISPIN        = 2
    MAGMOM       = 1.0 -1.0  1.0 -1.0  \
                1.0 -1.0  1.0 -1.0  \
                1.0 -1.0  1.0 -1.0  \
                1.0 -1.0  1.0 -1.0  \
                16*0.0
    LORBIT       = 11
    LMAXMIX      = 4 
    

This run provides the reference ground state charge density and wavefunction (CHGCAR, WAVECAR).

3. NSCF Step

Copy the following from the ground state directory:

• POSCAR
• POTCAR
• INCAR
• KPOINTS
• WAVECAR and CHGCAR

Then add the following line to INCAR for the NSCF step:

ICHARG = 11

Also append these LDA+U parameters (example for Ni with U and J = 0.10 eV, where only the first Ni species receives the correction):

    LDAU         = .TRUE.
    LDAUTYPE     = 3
    LDAUL        = 2 -1 -1           # Apply to d-orbitals of Ni 1
    LDAUU        = 0.10 0.00 0.00   # U parameter for Ni 1
    LDAUJ        = 0.10 0.00 0.00   # J parameter for Ni 1
    

In this example, the correction is applied to the first atom (Ni, 3d orbital). Adjust accordingly for your system.

4. SCF Step

From the NSCF run, copy:

• POSCAR
• POTCAR
• INCAR
• KPOINTS
• WAVECAR

Remove ICHARG = 11 from INCAR for the SCF run.

5. File Setup Note

The file setup described above is crucial for proper Hubbard U calculations. The POSCAR and POTCAR files must reflect the distinction between atomic species that receive different treatments, even if they are chemically identical atoms.

    TITEL  = PAW Ni 02Aug2007
    TITEL  = PAW Ni 02Aug2007
    TITEL  = PAW O  22Mar2012
    AFM  NiO
    4.03500000 
    2.0000000000   1.0000000000   1.0000000000 
    1.0000000000   2.0000000000   1.0000000000 
    1.0000000000   1.0000000000   2.0000000000 
    1 15   16
    

6. Parameter Sweep

Repeat NSCF and SCF runs for the following U = J values:

0.15, 0.20, 0.05, 0.00, -0.05, -0.10, -0.15, -0.20 eV

7. Analysis

After obtaining the number of localized electrons (e.g., f or d orbital) for each run, perform linear response fitting using Python:

    import numpy as np
    import matplotlib.pyplot as plt

    # Applied potential (eV)
    J = np.array([-0.1, -0.05, 0, 0.05, 0.1, 0.15, 0.2])

    # Number of f-electrons (example data)
    NSCF = np.array([2.182, 2.345, 2.517, 2.697, 2.885, 3.08, 3.282])
    SCF  = np.array([2.448, 2.461, 2.520, 2.533, 2.548, 2.563, 2.575])

    # Fit linear response (n vs J)
    s_n, intercept_n = np.polyfit(J, NSCF, 1)
    s_s, intercept_s = np.polyfit(J, SCF, 1)

    U = 1.0/s_s - 1.0/s_n

    # Generate fits
    J_fit = np.linspace(J.min(), J.max(), 200)
    fit_NSCF = s_n * J_fit + intercept_n
    fit_SCF  = s_s * J_fit + intercept_s

    # Print results
    print(f"chi0 (NSCF slope) = {s_n:.4f}")
    print(f"chi  (SCF slope)  = {s_s:.4f}")
    print(f"U (eV)            = {U:.4f}")

    # Plot (like VASP wiki)
    plt.figure(figsize=(7,5))
    plt.scatter(J, NSCF, color='black', label='NSCF', marker='o', facecolors='none', s=70)
    plt.scatter(J, SCF,  color='red', label='SCF',  marker='s', s=60)
    plt.plot(J_fit, fit_NSCF, color='black', linestyle='-', label=f'NSCF fit (χ₀={s_n:.3f})')
    plt.plot(J_fit, fit_SCF,  color='red', linestyle='-', label=f'SCF fit (χ={s_s:.3f})')
    plt.xlabel("Applied potential J (eV)")
    plt.ylabel("Number of f-electrons")
    plt.title("Linear Response for Hubbard U Determination of UGe₂")
    plt.grid(True, linestyle=':')
    plt.legend()
    plt.tight_layout()
    plt.show()
    

Example output:

    chi0 (NSCF slope) = 3.6700
    chi  (SCF slope)  = 0.4379
    U (eV)            = 2.0114
            

This method extracts the on-site Coulomb interaction U via linear response by comparing self-consistent and non-self-consistent changes in orbital occupations with applied potential perturbations.

Work function

Check back soon!

Spin-Orbit Coupled DFT+U (LSDA+U) Calculation

This tutorial demonstrates how to perform spin-orbit coupling (SOC) calculations with DFT+U corrections, particularly useful for heavy elements with strong spin-orbit interactions like actinides and lanthanides.

1. 1: Initial Collinear SCF Relaxation

First, perform a standard collinear spin-polarized structural relaxation to obtain the optimized geometry.

INCAR: [1] General Parameters + [2] Structure Optimization + [10] Spin Polarization + [13] GGA+U. Adjust MAGMOM and the LDAU lists to your species.

# INCAR for SOC step 1: [1] + [2] + [10] + [13]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [2] Structure Optimization
EDIFFG   = 1E-6        # convergence criterion for the ionic loop (pos = energy; neg = forces)
NSW      = 1000        # maximum number of ionic steps
IBRION   = 1           # ionic update: 1 = quasi-Newton, 2 = conjugate gradient (for bad guesses)
ISIF     = 3           # freedom: 2 = ions; 4 = ions + shape; 3 = ions + shape + volume
MAXMIX   = 80          # max steps stored in Broyden mixer (negative: resets after each ionic step)
ISMEAR   = 1           # smearing: -5 = insulators/semiconductors/DOS; 0 = large systems (SIGMA = 0.05); 1 = metals (SIGMA = 0.2)

## [10] Spin Polarization / Magnetism
ISPIN    = 2           # spin-polarized calculation
MAGMOM   = 4*0 4*1     # starting moment per atom, in POSCAR order (here 4 atoms at 0, then 4 at 1)

## [13] GGA+U
LDAU     = .TRUE.
LDAUTYPE = 1           # 1 = Liechtenstein et al.; 2 = Dudarev (only U - J is meaningful)
LDAUPRINT = 0
LDAUL    = -1 2 -1 -1  # orbital to correct, one per species in POSCAR order: -1 = none, 2 = d, 3 = f
LDAUU    = 0 3.71 0 0  # U (eV), one per species
LDAUJ    = 0 0 0 0     # J (eV), one per species
LASPH    = .TRUE.

2. 2: Collinear SCF Calculation

Run a collinear spin-polarized calculation (e.g., DOS calculation) to obtain CHGCAR and WAVECAR files using vasp_std. Use the following INCAR template (example for UTe₂):

# General Parameters/Electronic Relaxation
PREC = Accurate  # Precision
ENCUT = 500      # cutoff energy for the planewave basis set in eV
EDIFF  = 1E-8    # convergence condition for the electronic loop
NELM = 400       # maximum number of electronic iterations
NELMIN = 4       # minimum number of electronic iterations
KPAR = 1         # no. of nodes (max. no. of k-points)
NCORE = 4        # no. of cores per node
LWAVE = .TRUE.   # Write WAVECAR (needed for SOC step)
ISPIN = 2        # Spin-polarized calculation (required for SOC)
MAGMOM = 4*1 8*0 # Initial magnetic moments (adjust for your system)

# DFT+U parameters (example for Uranium f-electrons)
LDAU = .TRUE.
LDAUTYPE = 1     # Liechtenstein method
LDAUL = 3 -1     # 3 for U (f orbitals); -1 for Te (not applied)
LDAUU = 3.42 0   # U parameter in eV
LDAUJ = 0.51 0   # J parameter in eV
LMAXMIX = 6      # For f-electrons

# DOS calculation settings
NEDOS = 1001     # number of DOS gridpoints
LORBIT = 11      # DOSCAR and PROCAR are written
ISMEAR = 1       # smearing for metals
            

3. 3: Non-collinear SOC Calculation (not relaxing)

Perform the spin-orbit coupling calculation using vasp_ncl (non-collinear version of VASP). Copy the following files from Step 2:

• INCAR
• POSCAR
• POTCAR
• KPOINTS
• CHGCAR

Modify the INCAR file by adding the following SOC parameters:

# Add to existing INCAR from Step 2 and comment out the DOS-specific tags (not optimizing here), although I'm not 100% certain if they must be enabled off before calculating for the energy:
MAGMOM = 1 0 0 1 0 0 1 0 0 1 0 0 24*0  # x y z components of magnetic moments, adjust as per SAXIS
LSORBIT = .TRUE.  # Enable spin-orbit coupling
ISYM = 0          # Disable symmetry (recommended for SOC)
SAXIS = 1 0 0     # Quantization axis (easy axis)
ICHARG = 11       # Read charge density from collinear calculation
            

Important: Remove WAVECAR before running the SOC calculation, as the spin density is read from CHGCAR but the wavefunctions need to be recalculated for the non-collinear case.

4. 4: Probe Magnetic Anisotropy (optional)

To investigate magnetic anisotropy, run separate SOC calculations with different quantization axes (SAXIS) and see which is most favorable (lowest energy). Run three calculations with:

• SAXIS = 1 0 0 (a-axis, often the easy axis)
• SAXIS = 0 1 0 (b-axis)
• SAXIS = 0 0 1 (c-axis)

Compare the total energies to determine the magnetic easy axis and anisotropy energy.

5. 5: Running the Calculations

Execute the calculations using the appropriate VASP executable:

• Collinear calculations: mpirun -np [cores] vasp_std > output_collinear &
• SOC calculations: mpirun -np [cores] vasp_ncl > output_soc &
6. 6: Analysis

After completing the calculations:

• Compare total energies between different SAXIS orientations (optional)
• Analyze the magnetic moments from OUTCAR
• Check for changes in electronic structure via DOS analysis
• Evaluate the impact of SOC on the ground state magnetic configuration
• You can run DOS, band structure calculations with SOC tags enabled as above from step 3, just add necessary DOS/band structure tags as needed. Spin contributions are not included in the outputs.
7. Important Notes
• SOC calculations require significantly more computational resources than collinear calculations
• The MAGMOM tag in SOC calculations specifies x, y, z components of magnetic moments
• Always start from a converged collinear calculation
• For heavy elements, SOC can dramatically affect the electronic structure and magnetic properties
• Consider using hybrid functionals if standard DFT+U doesn't capture the physics correctly

Case study: Spin-orbit coupled DFT+U calculation of highly correlated uranium ferromagnetic superconductor system URhGe

This case study demonstrates the calculation of a highly correlated uranium ferromagnetic superconductor system URhGe using spin-orbit coupling combined with DFT+U methods. This is particularly important for actinide compounds where both strong electron correlation and spin-orbit coupling effects are significant.

Pseudopotentials

For this calculation, I used standard pseudopotentials

Part 1: Rh (d) U Correction

This part covers the calculation of Hubbard U correction for Rh d-electrons in the URhGe system. The tutorial video above demonstrates the methodology for determining the appropriate U value for Rh atoms.

Part 2: U (f) U Correction - Coming soon!

This section will cover the calculation of Hubbard U correction specifically for uranium f-electrons, which is crucial for accurately describing the strong electron correlation effects in actinide compounds.

Part 3: DFT+U Collinear Calculation - Coming soon!

Perform standard DFT+U calculations using the determined U values from Parts 1 and 2. This includes:

• Structural relaxation with DFT+U corrections
• Electronic structure calculations (DOS, band structure)
• Magnetic moment analysis

Part 4: DFT+U SOC Noncollinear Calculation - Coming soon!

Final calculations incorporating both DFT+U corrections and spin-orbit coupling effects:

• Noncollinear magnetic calculations with SOC
• Magnetic anisotropy analysis
• Final electronic structure with full relativistic effects

Note: For detailed steps on implementing Hubbard U correction, refer to the Hubbard U Correction section above. For spin-orbit coupling methodology, see the Spin-Orbit Coupled DFT+U section.

Fermi Surface

This section covers how to calculate and visualize Fermi surfaces using VASPKIT and PyProcar, which is essential for understanding electronic transport properties and band topology in metals.

1. 1: Initial Collinear SCF Relaxation

First, perform a standard collinear structural relaxation (can be spin-polarized) to obtain the optimized geometry: [1] + [2], optionally with [10] Spin Polarization and [13] GGA+U (see SOC).

2. 2: Collinear SCF Calculation

Run a collinear spin-polarized calculation (e.g., DOS calculation) to obtain CHGCAR and CHG (CHG may not be needed) files using vasp_std.

INCAR: the DOS setup ([1] + [3], [2] commented out), plus the three tags below. Comment out ISMEAR and LORBIT in [3] so they aren't set twice.

# INCAR for Fermi surface step 2: [1] + [3]
# (all other template sections commented out or deleted)

## [1] General Parameters
PREC     = Accurate    # precision
ENCUT    = 500         # plane-wave cutoff energy (eV)
EDIFF    = 1E-8        # convergence criterion for the electronic loop
NELM     = 150         # maximum number of electronic iterations
NELMIN   = 4           # minimum number of electronic iterations
KPAR     = 2           # k-points run in parallel (comment out for Bader)
NCORE    = 4           # cores working on each orbital
#LWAVE = .FALSE.   # skips writing WAVECAR; keep commented for DOS, band structure and COHP

## [3] DOS and Band Structure
NEDOS    = 1001        # number of DOS grid points
# if you add ICHARG = 11, a previous CHGCAR must be in the folder

## added to [1] for the Fermi surface
ISMEAR   = 0           # Gaussian smearing
ICHARG   = 11          # Read charge from CHGCAR
LORBIT   = 11          # Write orbital projections (for orbitals)

3. 3: Generate KPOINTS for Fermi Surface Calculation

Run VASPKIT to generate KPOINTS for the Fermi Surface calculation:

vaspkit
261 # Option for Fermi Surface KPOINTS generation
4. 4: Run the Calculation

Execute the VASP calculation with the generated KPOINTS file to obtain the electronic structure data needed for Fermi surface visualization.

5. 5: Export Required Files

After the calculation completes, export the following files for Fermi surface analysis:

• PROCAR - Contains orbital projections and band information
• KPOINTS - K-point mesh information
• POSCAR - Crystal structure (not needed for PROCAR but useful for reference)
• OUTCAR - Contains Fermi energy and other calculation details
6. 6: Visualization with PyProcar

Use PyProcar to visualize the Fermi surface. See the following PROCAR tutorials for different visualization options:

• 3D Fermi Surface: Getting Started with 3D Fermi Surfaces
• 2D Fermi Surface: Getting Started with 2D Fermi Surfaces
• 2D Fermi Surface with Spin: Dealing with Spin in 2D Fermi Surfaces
7. Important Notes
• Fermi surface calculations require dense k-point meshes for accurate results
• The quality of the Fermi surface depends heavily on the k-point sampling
• For spin-polarized systems, separate Fermi surfaces for spin-up and spin-down channels can be analyzed
• Ensure that the Fermi energy is correctly identified from the OUTCAR file

Phonons

Check back soon!

Processor Considerations

1. More processors doesn’t necessarily mean faster computations!
• Increasing the number of processors may not significantly accelerate the process due to several factors.
• Efficiency of parallelization parameters: It's important to optimize parameters like NCORE. For instance, setting NCORE=8 can enhance parallel performance when using 72 processors.
• Hardware limitations: The computer I utilized has 44 physical cores. Utilizing hyperthreading, the difference between running 2 threads versus one core is only about a 15% improvement.
• Diminishing returns: In this scenario, utilizing more than 32 cores may not substantially expedite the computation.
• Despite these considerations, it's worthwhile to experiment with different processor configurations to assess performance.

Troubleshooting errors

1. ERROR FEXCP: supplied exchange-correlation table is too small, maximal index: 5536
   • Solved by setting the following in [2] Structure Optimization (replacing the values there):
   • IBRION=1
   • ISIF=2
2. internal error in FOCK_FCC: number of k-points incorrect 5 18 2 2 2
   • Solved by ensuring KPOINTS are divisible by the grid reduction factors:
   • NKREDX = 2
   • NKREDY = 2
   • NKREDZ = 2
   • For example, a 5 x 5 x 5 mesh will not work. Use 4 x 4 x 4, etc.
   • Solution found here.

Further debugging resources

• Common VASP errors and solutions
• Materials Project Google Group: VASP errors
• VASP Forum: Error discussion 1
• VASP Forum: Error discussion 2
• VASP Forum: Error discussion 3

Calculating NBANDS (for metals)

• Set this in [5] COHP Analysis of the INCAR.
  • # of bands = (# of electrons/2) + (# of ions*2)
  • Running grep NELECT OUTCAR after an optimization will give you the number of electrons, and grep NIONS OUTCAR will give you the number of ions
  • Running grep VRHFIN POTCAR may also give you the number of bands (not certain, have to double check)