Using VASP for density functional theory (DFT)

Handy Terminal Commands

Login: ssh [username]@[serverIP]

Password: [password]

Copying a file from your local machine to the destination server: scp ~/[local file path] [username]@[serverIP]:~

• Note that moving files between your local computer and the server must be done within a terminal without being logged into the server

Copying a file from the server to your local machine: scp [username]@[serverIP]:~/[server file path] ~/[local file path]

Copying a folder from your local machine to the destination server: scp -r ~/[local file path] [username]@[serverIP]:~

Check CPU% of server and status of files: top or htop

Finding a job directory from PID #: pwdx [PID#]

See files in the working directory: ls

Delete files: rm [filename]

Delete folders + contents: rm -r [foldername]

Make a directory: mkdir [foldername]

Copy folders + contents: cp [foldername] [destinationfolder]

View & edit a file: vim [filename]

Rename a file: mv [oldfilename] [newfilename]

• To enter insert mode, click esc
• To save and exit, type esc (to leave insert mode) + :wq + enter
• To exit without saving after typing, type esc (to leave insert mode) + :q! + enter

Before You Start: Creating the Necessary Files

POSCAR file

1. Open CIF files with VESTA, which you can download here.
2. Export structure as fractional coordinate (rename the files as POSCAR in a folder with whatever naming convention you prefer)
   • On Macbook, it can only saved as POSCAR.vasp so you can't just delete it. You can only remove the .vasp with the command mv POSCAR.vasp POSCAR
3. Move POSCAR files from local computer to your folder on the destination server with terminal/powershell using the above keys
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.

POTCAR file

1. Check order of POSCAR molecules for the potcreate command
   • For example: potcreate Si Fe Ce in command line, which matches the order of molecules in the POSCAR file
2. Check available element potentials like this: potshow Si
3. 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

INCAR file

1. Copy the INCAR from this template below, or download it here.
### INCAR-Template for VASP Calculations - Martin Changed 2015-12-22
## 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 = 150 # maximum number of electronic iterations
NELMIN = 4 # minimum number of electronic iterations
KPAR = 2 # no. of nodes (max. no. of k-points)
NCORE = 4 # no. of cores per node
#LWAVE = .FALSE.# Do not write WAVECAR
## Structure Optimization (or static for NSW=0)
EDIFFG = 1E-6 # convergence criterion for ionic loops (pos=energy; neg=forces)
NSW = 1000 # number of ionic steps
IBRION = 1 # ionic update: 1:QN, 2:CG (for bad guesses)
ISIF = 3 # freedom: 2=ions; 4=ions,shape; 3=ions,shape,vol
MAXMIX = 80 # max number steps stored in Broyden mixer, when negative: resets after ionic step
ISMEAR = 1 # smearing: -5-insul/scond/DOS 0-large systems (with SIGMA 0.05) 1-metals (with SIGMA 0.2)
## DOS and Band Structure Calcuations (for ICHARG=11 prev. CHGCAR required)
#NEDOS = 1001 # number of DOS gridpoints
#LORBIT = 10 # DOSCAR and PROCAR are written
#ISMEAR = -5 # smearing: -5-insul/scond/DOS 0-large systems 1-metals (with SIGMA)
## Bader Charge and Electron localization function
#LAECHG = .TRUE.
#LELF=.TRUE.
#NPAR = 1
## COHP Analysis (via LOBSTER)
#LWAVE = .TRUE.
#NSW = 0
#ISYM = -1
#NBANDS = 600 ! should equal the total number of orbitals
#LORBIT = 12
#ISTART = 0
## Phonon Dispersion via finite displacements
#EDIFF = 1E-8 # delete other EDIFF tags!!!
#ISMEAR = 0
#SIGMA = 0.01
#ADDGRID = .TRUE.
## Hybrid Calculations
LHFCALC=.TRUE. # use of HF/DFT hybrid functionals
AEXX = 0.25 # fraction of exact exchange
AGGAX = 0.75 # fraction of GGA
ALGO = Damped # electronic minimisation algorithm (damped velocity friction algo.)
TIME = 0.05 # time step for ALGO
PRECFOCK = FAST # controls FFT grid for the exact exchange
HFSCREEN = 0.2 # determines the range separation parameter
NKREDX = 2 # grid reduction factors
NKREDY = 2 # grid reduction factors
NKREDZ = 2 # grid reduction factors
## Stress Tensor Calculation
#ENCUT = 600 # delete other ENCUT tags!!!
#IBRION = 6
#ISIF = 6
#NSW = 1
#NFREE = 2
## additionally for DFPT:
#IBRION = 8
#NWRITE = 3
#NSW = 1
## Supercell (band offset) calculations
#LVHAR = .TRUE.
#AMIN = 0.1
## Born Effective Charge
#LEPSILON = .TRUE. ! dielectric matrix
#LOPTICS = .TRUE. ! DFPT including empty bands for hybrid functionals
## Spin-orbing Coupling
#LSORBIT = .TRUE.
#LNONCOLLINEAR = .TRUE.
##GGA+U, Coulomb interaction energy
#LDAU = .TRUE.
#LDAUTYPE = 1 !1 = Liechtenstein et al.; 2 = Dudarev's approach, only the difference (U-J) is meaningfull
#LDAUPRINT = 0
#LDAUU = 0 3.71 0 0
#LDAUJ = 0 0 0 0
#LASPH = .TRUE.

• Adjust accordingly per your calculation of course. The above template indicates a structure optimization, as it's not commented out with '#'

KPOINTS file

1. Generate the file using VASPKIT, as seen in the video here .
2. You can also use VASPKIT by typing vaspkit in the terminal, then 01 for VASP Input-Files Generator, and 102 for KPOINTS file.
3. This file is an attempt to solve Schrodinger's equation… which is technically an infinite attempt as you have a guess and set up the cutoff and wanted precision in POSCAR file
   • K-Points: This line indicates the type of k-points generation scheme being used. In this case, it's set to "Gamma", meaning a single k-point at the center of the Brillouin zone.
   • 0: This value typically denotes the number of k-points to be generated. Here, "0" implies a single k-point.
   • Gamma: Specifies the center of the Brillouin zone where the k-point is placed.
   • 8 2 8: Defines the mesh size, i.e., the number of k-points along each reciprocal lattice vector direction. "8 2 8" indicates 8 k-points along the first lattice vector, 2 along the second, and 8 along the third.
   • 0 0 0: Specifies any fractional shifts applied to the k-point grid. In this example, there are no shifts applied.

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

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.
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)
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, indicated by EDIFF in 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

1. If a certain potential cannot converge no matter what, you need to enable the section for the hybrid calculations in the INCAR file. The optimization will take much longer, up to weeks in some cases.
2. In case there’s a problem with convergence, the way to fix that would be changing TIME = 0.2 to 0.1 in the INCAR file.
3. Setting an easier convergence criterion, namely 10^-6 and 10^-5 for electronic and ionic loop respectively may help with stubborn structures.

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

1. Uncomment DOS & band structure and General Parameters (never LWAVE = .FALSE.) in INCAR as we need the WAVECAR file for this.
2. Non-spin polarized DOSCAR
3. Spin-Polarized and Non-Spin-Polarized :
   • ISPIN=1 (non-spin-polarized) and ISPIN=2 (spin-polarized) calculations.
4. 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).
5. Run with mpirun -np 16 vasp_std>Dos1& (Oliynyk lab - use sbatch job after setting job file)
6. Copy DOSCAR files to local desktop for analysis when converged
7. 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)
   • Add the following to the INCAR:
   • ISPIN = 2
   • MAGMOM = ie.) 4*0 4*1 OR 0 0 0 0 1 1 1 1
   • This can only be full integers
   • 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
   • From there, plot the same elemental contributions in wxDragon like before

Electron Localization Function

1. Uncomment Bader Charge and Electron localization function and General Parameters in INCAR
2. Run with mpirun -np 16 vasp_std>Elf& (Oliynyk lab - use sbatch job after setting job file)
3. Import ELFCAR file into local computer and drag into VESPA to view

Bader Charges

1. Uncomment Bader Charge and Electron localization function and General Parameters in INCAR, comment out KPAR
2. Run with -np 8 vasp_std>Bader& (Oliynyk lab - use sbatch job after setting job file)
3. When converged run chgsum.pl AECCAR0 AECCAR2
   • This combines core and valence electrons
4. type bader CHGCAR -ref CHGCAR_sum
   • Looks like this:
   
5. Import ACF.dat to local machine or just view it in the terminal with vim command
6. Run grep ZVAL POTCAR
   • Looks like this:
   
7. Run vim POSCAR if you forget the order of the molecules to determine which charge is which molecule
8. Use above values as the initial charge, see the VASP output in ACF.dat to determine difference, or actual charge

Crystal Orbital Hamiltonian/Overlap Population

1. Uncomment COHP Analysis (via LOBSTER) and General Parameters in INCAR
2. 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

3. type mpirun -np 20 vasp_std>COHP& (Oliynyk lab - use sbatch job after setting job file)
4. type lobster
5. Note: If you do not set NBANDS 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.
6. Import COOPCAR.lobster, COHPCAR.lobster, to local machine as needed for plotting
7. To plot COHP/COOP, save the files (COOPCAR.lobster, COHPCAR.lobster, etc.) in a .zip archive and upload to COHP-COOP-GUI.
8. See GitHub for more info and settings.

Band Structure

1. Uncomment DOS & band structure and General Parameters (never LWAVE = .FALSE.) in INCAR as we need the WAVECAR file for this.
2. Then 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. Use the standard INCAR settings for structure optimization with DFT+U parameters.

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. Use the standard INCAR settings for structure optimization (Can use DFT+U as well, 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. Add the following to the INCAR's general parameters:

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 adding the following to the INCAR:
   • 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