Disclaimer

Error Handling

Input Flags

Geometry Optimizations

• Don’t waste your time using super high-accuracy settings on a structure far from the local minimum. Do an initial optimization with “looser” settings (e.g. in terms of ENCUT and k-point grid) and then refine it with the desired parameters.
• When performing geometry optimizations that involve changes in the cell shape and/or volume, always ensure that ENCUT >= 1.3*ENMAX to prevent  Pulay stresses . For structures containing carbon, this often means using a minimum cutoff of 520 eV.
• The first step in a geometry optimization will generally have the highest number of SCF iterations. It is okay if that first step does not converge electronically within the limits of NELM. In fact, it is better to have the first step reach NELM instead of running for many hundreds of SCF iterations.
• The general rule-of-thumb is to set EDIFFG to 0.05 eV/Å or lower for optimizations (note: that means EDIFFG=-0.05 since negative signs are used for force tolerances). I personally recommend 0.03 eV/Å or lower. For large, flexible materials, any value greater than 0.03 eV/Å is likely too high. Remember that in VASP, the EDIFFG flag must be a negative value if you want a convergence threshold based on the maximum net force, not the change in total energy.
• When performing a full optimization of atomic positions and cell shape/volume, it is almost always best to do this in stages. It is often wise to start with a relaxation of the atomic positions (ISIF=2) followed by a full volume relaxation (ISIF=3). This will significantly reduce the chance of running into convergence issues.
• If your material can only be modeled with vacuum space as part of the simulation unit cell (e.g. a surface slab), never use ISIF=3 to optimize the cell volume, as it will simply reduce/eliminate the vacuum. You need to vary the lattice parameters manually and find the global energy minimum that way.
• When performing optimizations, you are probably better off setting ISYM=0 unless you are certain you wish to constrain the symmetry. Even though this will increase the computational cost by not using available symmetry, it allows the system to reach local minima that would not otherwise be accessible. In addition, it prevents you from having to worry about errors related to SYMPREC. That being said, some symmetry issues sometimes still arise despite having ISYM=0. To resolve these, set SYMPREC=1.0e-8.
• If you wish to stop a job but want it to still output any WAVECAR or CHGCAR files, create a file in the working directory named STOPCAR and have a line that reads LSTOP=.TRUE.(this will cause the job to stop on the next ionic step, and all restart files will be written)

Choice of Geometry Optimization Algorithms

• In choosing an optimization algorithm for finding local minima, I generally recommend starting with the conjugate gradient (CG) algorithm (IBRION=2) because it is very robust and you do not have to worry about tweaking POTIM.
• In large, flexible materials with many degrees of freedom, the CG optimization algorithm (IBRION=2) oftentimes results in a bracketing error once it gets relatively close to the local minimum (search for ZBRENT: fatal error in bracketing in the standard output file). This occurs because the potential energy surface is very flat, and the CG algorithm implemented in VASP is based on the energy differences. One option to fix this is to copy the CONTCAR to the POSCAR and tighten EDIFF to 1e-6 , but a more reliable option is to use a force-based optimizer. Of these, I’d recommend FIRE as implemented with VTST (IBRION=3, IOPT=7). I have found that FIRE is generally more robust than the QN (IBRION=1) method.
• If the starting structure has extremely high forces on some atoms, make sure you use an appropriate optimizer, at least for the first few steps. I strongly recommend a force-based optimization algorithm with a robust line-search algorithm. For cases where the forces are so high that the structure “explodes” within a few iterations, I recommend using ASE’s BFGSLineSearch algorithm until max|F| < 10 eV/A or so. Then continue with your favorite optimizer.
• Generally, for NEB and CI-NEB calculations,  the L-BFGS algorithm (IOPT=1) implemented in VTST is the fastest. For the dimer method, the force-based CG method in VTST (IOPT=2) is recommended. However, if you are having trouble in either case, I suggest switching to the FIRE algorithm (IOPT=7) with the default settings. It is a bit slower, but it is especially useful in troublesome cases of convergence.

Electronic Energy Convergence

• For insulating materials, SCF convergence is greatly accelerated by using ALGO=All.This has the added benefit that you don’t have to worry about any of the mixing tags. This algorithm is also recommended when using the M06-L meta-GGA functional, regardless of system type.
• I recommend setting LWAVE=.TRUE. for most jobs, unless you are using an ASE-based optimizer or running many short jobs in succession. There is a large file I/O cost, but the WAVECAR serves as an excellent starting guess for continuation jobs, if needed, and can be deleted when you’re done. I don’t necessarily recommend setting LCHARG=.TRUE. for optimizations. The initial charge density can be calculated from the previous WAVECAR. 
• Be careful about obtaining final energies from calculations using real-space projectors (e.g. LREAL=Auto) even for large systems. This can impact energies somewhat (even energy differences) depending on what you wish to study. If VASP recommends in the OUTCAR file that your system is large and that you will save time by using real-space projectors, I strongly recommend running your calculations with LREAL=Auto and then doing a final optimization with LREAL=.False. The optimization with the reciprocal-space projectors will only require an additional 1-5 steps on average, so you will still get the speedup of using real-space projectors.
• While the VASP manual suggests not setting the maximum number of SCF iterations (NELM) above 60, there are many materials (e.g. MOFs) where a higher value is needed for the first few steps. I generally set NELM=150 when studying very large materials.
• For materials that have charge sloshing or difficult convergence issues with the SCF, I recommend setting NELMIN to 4-6 so that accurate energies and forces are obtained.
• If you do not explicitly set ISTART or ICHARG, VASP will detect if a WAVECAR is present and use this for the initial wavefunction and charge density if it exists. This is arguably better than explicitly setting ISTART in your INCAR because it prevents VASP from crashing if a WAVECAR is not present (e.g. because the previous job crashed).
• For insulating materials (or materials that have unknown band gaps), I recommend using Gaussian smearing (ISMEAR=0) since it is appropriate for both conductors and insulators, although there are better choices for conducting materials.
• For insulating materials, I recommend settings SIGMA to 0.01 eV as a reasonable initial test. Always check to make sure that extrapolation back to T = 0 K from the finite-temperature approximation is close to the fictitious free energy, as listed in the OUTCAR file.
• It is generally good practice to set PREC=Accurate. This provides reasonable estimates for the integration grid (NGX, NGY, NGZ). It is often stated the meta-GGA functionals require tight integrations grids to achieve appreciable convergence. I have found that PREC=Accurate is sufficient even for these troublesome functionals.
• If you are going to write out or read in a charge density, always set LMAXMIX=4 if you have d-block elements or LMAXMIX=6 if you have f-block elements. While the effect will be negligible if you use ICHARG=1, if you run a non-self-consistent calculation with ICHARG=11, omitting LMAXMIX will cause the total energy and computed properties to differ. This is especially true for spin-polarized systems since the magnetic moments are not fully reproduced.

Spin Polarization

• When studying systems with magnetic moments, set LORBIT=11 so that you can view the converged magnetic moments for each atom.
• If you are continuing a job from a previous WAVECAR, you do not need to set the MAGMOMS in the INCAR file so long as ISPIN=2 (you can if you’d like, but VASP will ignore them).
• If you are unsure what to set as the magnetic moment for a metal atom, consider setting the MAGMOM value to the number of anticipated unpaired electrons. If this too is completely unknown, setting a value of 5 for d-block and 7 for f-block elements is typically okay as a first guess. If this approach is taken, it would be wise to use the converged structure to also test other spin initializations. I would also test an initial MAGMOM of 0.1 to see if this converges to a system with no spin, back to the same magnetic moment as the high-spin initialization, or something else entirely.

Spin/Charge Densities

• Never generate AECCAR files (e.g. for a Bader analysis or DDEC charges) during a geometry optimization. VASP writes the AECCAR0 file for the input geometry but the AECCAR2 file for the converged geometry, so they are not compatible. Always do a separate single-point energy calculation (nsw=0) when setting LAECHG=.True.

Pseudopotentials

• Always consider using the VASP-recommended PAW potentials, shown  here . Note that Li_sv has an ENMAX of 499 eV in the 5.4 version of the PAW_PBE pseudopotentials, but all the rest should have ENMAX < 400 eV. There is now an option in ASE to automatically choose the VASP-recommended default pseudopotentials. Just set setups='recommended'.
• For DFT+U, hybrid functionals, and meta-GGA calculations, it is  recommended  to set LASPH=.TRUE. If you use the LASPH=.TRUE. flag, make sure you include it for gas-phase species as well! It should be considered as part of your model chemistry.
• A subtle point is that, because the non-spherical contributions to the gradient corrections inside the PAW spheres will change with ENCUT, it will be hard (if not impossible) to converge an absolute energy with respect to the plane-wave kinetic energy cutoff if LASPH=.TRUE. Instead, one should look at convergence of relative energies or do a convergence test with LASPH=.FALSE.

Parallel Performance

• The rule-of-thumb is that the number of processors you should use scales reasonably well with the number of atoms in the system. This holds fairly well for metallic systems where there are many electron-rich elements, but can be a slight over-estimate for organic or organometallic systems where there are many C and H atoms.
• Whenever possible, use the gamma-point only VASP executable, as it runs up to 1.5 times faster than the standard executable.
• I recommend using NCORE in place of NPAR since it automatically adjusts based on the number of nodes. The optimal value for NCORE strongly varies based on the computing environment. A good first-guess on many compute clusters is to try setting NCORE to the number of processors in a given node. Regardless, always make sure that the number of processors per node is divisible by NCORE.

Third-Party Utilities

VTST

• If upon restarting a dimer calculation you find that the torque and angles are higher than where they left off, make sure that you are using a VASP build with VTST 3.2 or newer. See  here .
• The usage instructions for the VTST dimmins.pl script are incorrect. It should be dimmins.pl POSCAR MODECAR displacement, where the POSCAR and MODECAR are the resulting files from doing  vfin.pl 
• The usage instructions for the VTST neb2dim.pl script are unclear. It should state that you must first run vfin.pl, copy the exts.dat file to the parent directory (where the new POSCAR files are written), and run neb2dim.pl from that parent directory.

Atomic Simulation Environment

• Whenever possible, refrain from using the ASE optimizers (or any external optimizer) with VASP calculations. If necessary, consider using the VaspInteractive calculator instead of the Vasp calculator, as this will reduce CPU time associated with starting and stopping VASP every ionic step. However, even if the VaspInteractive calculator is used in place of the Vasp calculator, VASP will generally require more electronic steps when run with an ASE optimizer than with an internal VASP or VTST optimizer. This is because VASP keeps an internal history of previous ionic steps in its mixer history and because VASP extrapolates the wavefunction and charge density from the previous step to the new positions. There is also the I/O overhead associated with reading/writing any restart files every ionic step.
• The above holds true for vibrational frequency calculations as well. That being said, it sometimes may be worth it to use ASE’s vibrations module even if there is added CPU overhead. The ASE-generated vibrations can be easily visualized using the ASE GUI, whereas VASP-generated vibrations are more difficult to visualize. The ASE vibrations module also saves restart files (.pckl files) for every displacement, so if the job crashes or exceeds the specified walltime, it is easy to continue where you left off.