2.10. MedeA VASP 5.4

---

download:pdf

---

2.10.1. Introduction and Capabilities:

Note

VASP [1] is applicable to bulk solids, surfaces, interfaces, molecules on surfaces, and molecules. VASP is a fast and highly reliable electronic structure method based on density functional theory (DFT) [2]. Together with the all-electron projector augmented wave potentials [3] VASP has the generality and accuracy of an all-electron method while maintaining the speed and advantages of a plane-wave method. Advanced simulations beyond DFT enable highly accurate predictions much beyond standard DFT precision.

MedeA’s graphical VASP user interface gives easy access to the relevant VASP parameters in a structured arrangement of topic panels. In addition, the MedeA interface provides a comprehensive and exhaustive set of defaults for all relevant parameters needed to run standard VASP calculations. For special settings, direct additions to VASP input files through MedeA are provided.

VASP provides basic properties such as total energies, optimized geometries, trajectories, band structure and density of states plots, charge density, potential, and magnetization data, charge analysis, optical spectra, zone center phonon frequencies, response tensors (dielectric, piezoelectric, and Born effective charge tensor), electric field gradients, hyperfine parameters, and NMR chemical shifts. In addition, MedeA modules use VASP as a DFT solver for the electronic total energy, interatomic forces and the stress tensor to give access to more complex properties derived from “simple” single point runs or structure optimizations. Examples are elastic constants, phonon spectra, thermodynamic functions, Fermi surfaces, electronic transport properties, effective masses, and transition states.

VASP 5.4 Fundamental Capabilities

Single Point

Total energy for a given fixed arrangement of atoms in a unit cell (atoms, molecules, surfaces, solids); Energy of formation, reaction, adsorption, or insertion; Relative phase stability at 0K.

Structure Optimization

Structure of bulk solids, surfaces, interfaces, molecules on surfaces by minimizing forces and energy, relaxing lattice parameters and/or internal degrees of freedom in the process: crystal structure determination, surface and interface relaxation, bond lengths, adsorption geometries, defects and vacancies

Ab Initio Molecular Dynamics

Evolution of a system of ions by calculating electronic structure and inter atomic forces from ab-initio and applying Newton mechanics for the ionic movement; Temperature effects through kinetic energy term; Structure determination for complex systems; Simulated annealing

Electronic Structure

Electronic band structures and density of states, charge densities, electronic localization function, total potential, magnetization densities, total charge density, total valence charge density, Bader charge analysis, work function, magnetic moments, bonding, optical spectra, electronics, dielectric tensor, piezoelectric tensor, Born effective charges, zone center phonon frequencies, electric field gradients and quadrupole coupling constants, hyperfine parameters, NMR chemical shifts

History

Note

VASP [4] is short for “Vienna Ab initio Simulation Package”. The program has originally been developed in the group of Prof. Jurgen Hafner, who was heading the Institute for Materials Physics at the University of Vienna, Austria. In the early days, Georg Kresse and Jurgen Furthmuller have been the key authors of VASP. Nowadays, the VASP code is continuously developed by Georg Kresse and Martijn Marsman and a dedicated team of co-workers of the VASP Software GmbH. The Wiki based VASP Manual provides an excellent description of the algorithms underlying VASP and a detailed compendium of computational options.

VASP can be considered as the culmination of many decades of worldwide efforts in electronic structure theory. Through the implementation [5] of the projector augmented wave (PAW) method [6] VASP combines the speed and elegance of plane-wave methods with key features of frozen-core all-electron methods.

The integration of VASP and MedeA is far beyond a simple graphical user interface: For example, the automatic calculation of elastic constants relies on years of experience usually gained during a Ph.D. length of time. Sound knowledge of group theory is required to create the specific distortions and combine the results into an elastic matrix. On a technical side the various competing options of potentials, integration methods, and cut off energies have to be rigorously tested before using them in industrial R&D projects.

MedeA’s graphical user interface is on top of a set of tested parameters, optimized for specific types of computation with very different levels of required accuracy. If needed, this interface gives experts access to the less frequently used features of VASP through viewing and editing capabilities of the standard VASP input files. Non experts find good default values of computational parameters and context-sensitive help.

A detailed, technical description of the underlying algorithms is found at https://www.vasp.at/wiki/index.php/The_VASP_Manual

Note on the automatic computation of the heat of formation for compounds:

MedeA VASP 5.4 does not compute automatically the heat of formation for a given compound, as MedeA VASP 5.2 and MedeA VASP 4.6 do, only VASP total energies are available. The vast number of different combinations of options to evaluate the interaction makes the approach pursued in the MedeA VASP 5.2 and 4.6 GUIs impracticable. The capability may be replaced by different procedures in the future.

The heat of formation of a compound is defined as the difference in enthalpy at room temperature between the compound and its constituent elements in their standard state. Contributions to this property are:

• The electronic term (difference of VASP total energies), temperature not considered
• The temperature dependent electronic term (generally very small)
• The zero-point energy term
• The contribution of lattice vibrations at room temperature

A first (often quite good) approximation is given by the electronic contribution.

2.10.2. The MedeA-VASP 5.4 Interface

Getting started

From the MedeA toolbar, select Tools >> VASP 5.4. A new menu entry, VASP 5.4 appears in the MedeA menu bar and remains there for the rest of the current MedeA session.

Load or select the system, for which you want to perform a VASP calculation. Note that the File >> Open structure from disk command allows the import of structural information from POSCAR/CONTCAR files. If these structure files were created by VASP 4.6 the corresponding POTCAR file must exist in the same directory, if created by VASP 5.2 or VASP 5.4, the POTCAR file is not required anymore. The structure information for VASP 4.6 is spread over two files: atomic coordinates in POSCAR/CONTCAR and element information in POTCAR. For VASP 5.2 and VASP 5.4 the element information is contained in POSCAR/CONTCAR as well. This is useful to import earlier VASP calculations (e.g. those run outside MedeA or from those of interrupted tasks). In addition, several other external file formats can be read: .cif, .car, .xtl, .xyz, .mol, .log

Select VASP 5.4 >> Run from the MedeA toolbar to bring up the VASP graphical user interface. In the VASP user interface, right-click on a text field to get context-sensitive help or to reset a parameter to its default value.

Interface panel overview

The MedeA VASP graphical user interface consists of a stack of panels or cards grouping the main parameters relevant for the setup of VASP runs.

For all input parameters, MedeA provides defaults that were chosen to yield acceptable precision while limiting the computational effort. It is strongly recommended to study convergence for each specific case before concluding any scientific results. However, in many cases the defaults provided give an excellent starting point for the structural and electronic properties of a given system. Converged total energies and heats of formation usually require a higher level of accuracy (so-called Standard 500 settings).

At the bottom of the MedeA VASP GUI there are five buttons to execute tasks affecting and taking into account all panels:

Run: Submits and runs all specified calculations with the chosen settings on the active JobServer and attached TaskServers

Close: Closes the graphical user interface and memorizes all settings for later use

Write input files: Writes representative VASP input files to a selected destination folder

Restore defaults: Restores MedeA VASP default values for all parameters and settings

Restore from job: Restores all parameters and settings from a previous job, to be selected from a browser panel

2.10.3. VASP Output Files:

During VASP runs, data is written to the task directory on a TaskServer machine. Upon successful completion of a computational job, all VASP tasks associated with the job are transferred back to the JobServer. You can access output files on the JobServer by browsing to the job directory where all data from completed runs are stored (in the MedeA menu click Jobs >> View and Control Jobs).

The most relevant output files are

• The Job.out summarizes the entire computational job from input parameters to the main results.
• The VASP.out or OSZICAR.out summarizes convergence information for geometry steps and electronic iterations of a given VASP task.
• The OUTCAR.out: is the detailed output from a given VASP tasks.

Note

For any given crystal symmetry, VASP performs the calculations on the primitive cell. If you are in doubt if the actual cell shape used by VASP corresponds to the one displayed by MedeA, check if a primitive cell exists (Symmetry panel) or simply choose to reduce symmetry to P1 in MedeA.

Example:

	Fm-3m	P1	P1 primitive
Cell displayed by MedeA
Cell used by VASP

2.10.4. Interface Description

Click on the tabulators Calculation, Functional/Potential, SCF, DOS/Optic/Tensors, Band structure, Advanced/Restart, etc. to go to the respective panel, which are described in detail in the following.

2.10.5. The Calculation Panel

Here you set the type of calculation to perform using VASP, you select the properties to be evaluated, you may impose external conditions such as solvation and pressure and possibly a charge state, you chose how to handle the interactions between atoms, and you specify the general setup and accuracy requirements for the calculations. To this end the panel has several submenus/fields:

• Type of calculation
• Properties
• Solvation, pressure, charge state, external electrostatic field
• Interaction
• General Setup

2.10.5.1. Type of Calculation

The MedeA VASP interface distinguishes 7 major types of calculations, namely Single point, Structure Optimization, Molecular Dynamics, Time-dependent hybrid / DFT, Quasiparticle Spectra (GW), Accurate Energy (ACFDT-RPA) and MT -Elastic Properties. The MT entry is visible only with a valid MT license.

Single Point: Performs an electronic structure calculation for the input geometry without relaxing any structural parameter.

Structure Optimization: Relaxes the atomic position and/or the cell parameters with or without constraints. Perform a full structure optimization to determine the bulk equilibrium structure at T=0K

Molecular Dynamics: Performs ab-initio molecular dynamics for dynamic properties or equilibrium states

Time-dependent hybrid / DFT: Evaluates accurate optical spectra including excitonic effects based on DFT or hybrid functional calculations for the electronic structure

Quasiparticle Spectra (GW): Calculates quasiparticle energies (excitation energies) employing the so-called GW approach. In addition, excitonic effects in the optical spectra can be obtained, based on the quasiparticle spectra by solving the Bethe-Salpeter equation.

Accurate Energy (ACFDT-RPA): Evaluates very accurately the total energy of the input structure (without geometry optimization) as a sum of the exact exchange energy and the correlation energy within the random phase approximation (RPA) employing the adiabatic connection fluctuation dissipation theorem (ACFDT).

MT - Elastic Properties: Computes elastic constants and other mechanical and thermodynamic properties (based on the Debye model)

2.10.5.1.1. Single Point

No further parameters required in this context field. However, you may want to check additional settings, in particular, you should know about the plane wave cutoff and the setting of the k-mesh (see SCF -Tab).

2.10.5.1.2. Structure Optimization

Selecting Structure Optimization activates the field for Structure Optimization parameters.

Available options are listed below, followed by recommended settings for standard tasks.

Relax atom positions: Atoms are moved until forces are smaller than the value in Convergence in eV/Ang

Allow cell volume to change: Varies volume while keeping constant ratios of a:b:c with unchanged cell angles

Allow cell shape to change: Varies the ratio a:b:c and changes the cell angles

Goal	Recommended Structure Optimization settings
Determine a bulk equilibrium crystal structure	Relax atom positions Allow cell volume to change Allow cell shape to change
Relax a surface (no in-plane relaxation); Find the equilibrium geometry of a molecule; Adsorb a molecule on a surface; Locally relax a structure around a vacancy/defect	Relax atom positions Allow cell volume to change Allow cell shape to change
Optimize a system under pressure; Allow in-plane adjustment during a surface calculation	Relax atom positions Allow cell volume to change Allow cell shape to change

Any combination of these three parameters can be chosen with one exception, which is not implemented in VASP: Unit cell volume and atom positions can only be simultaneously relaxed if also the unit cell shape is relaxed.

Update algorithm:

Conjugate Gradient (CG) is the default for standard structure relaxations minimizing forces and energy to find a local minimum of the total energy surface

RMM-DIIS is a Newton-Raphson based algorithm that converges faster than CG, if (and only if) the initial system is close to an extremum of the total energy surface. The algorithm is based on forces only, disregarding energies. Consequently, RMM-DIIS is able to converge to a saddle point or to a minimum, if the starting configuration is close enough to these extrema.

Convergence: Set an upper limit for the largest allowed residual force between any of the atoms in the unit cell. A value of 0.02 eV/\({\mathring{\mathrm{A}}}\) is reasonable for most calculations.

• High precision calculations may require 0.01 eV/\({\mathring{\mathrm{A}}}\) or even smaller residual forces.
• When reducing the criterion for the force convergence, you must use a lower value for the SCF convergence (see SCF panel), too: Try 1.0e-06 to 1.0e-08.

Maximum number of steps: Sets the maximum number of geometry steps to be executed before stopping. Roughly, the number of geometry steps can be of the same order of magnitude as the number of degrees of freedom present in the system. If the number of degrees of freedom is very large, one may consider using molecular dynamics and simulated annealing to find the minimum structures.

Trajectory file frequency: sets the number of animation frames to be written to disk during a geometry optimization. Default is 1 frame/geometry step. If set to 0 a trajectory file is not created.

Note: At present, the combination of only Relax atom positions and Allow cell volume to change checked are not handled by VASP.

2.10.5.1.3. Molecular Dynamics

In a molecular dynamics run the forces calculated in a given geometry step are used to update the atomic positions. The system dynamics, i.e. the ionic movements are subject to Newton mechanics while the forces acting on the ions are calculated from ab-initio using a self-consistent electronic density (Hellmann-Feynman forces).

Note

Note that the natural time step of an ab-initio molecular dynamics run is rather short compared to the time span required for a chemical reaction: For ab-initio dynamics, the typical time range is in the picosecond range! Thus results from a picosecond range ab-initio dynamics run need to be interpreted with care: If you observe a specific event in the analysis of an ab-initio molecular dynamics run, the event is most likely to have physical significance. However, the absence of events, for example, the absence of a diffusion jump from a molecular dynamics run cannot immediately be interpreted as a factual result. It may simply mean that the statistical sampling was too short.

Four types of molecular dynamics run are currently implemented in VASP and are accessible from the MedeA interface by the Ensemble menu. Their purpose and related settings are explained below:

Micro canonical (nVE): Molecular dynamics at a constant number of particles, n, constant volume, V, and constant free energy, E. The free energy consists of the electronic energy and the Madelung energy and kinetic energy of the ions. The parameters below apply for all ensembles:

Simulation time: The overall simulation time of the molecular dynamics run in femto seconds (fs)

Time step: Default is 4 femto seconds (fs), set time step to 1 fs if hydrogen is present

Temperature initial: Initial velocities are randomly attributed to all atoms according to a Maxwell-Boltzmann statistics at this chosen temperature

Temperature scaling (nVE): Simulated annealing to find energy minima for a complex structure with many degrees of freedom.

Temperature end: Final temperature at the end of the molecular dynamics simulation

Trajectory file frequency: For all ensembles this sets how frequent pair correlation functions and density of states are reported in the output file PCDAT. As a default this is done each geometry step. In simulated annealing mode this parameter specifies how frequently the kinetic energy (and thereby the temperature) is scaled to enforce the temperature gradient. If set to 0 a trajectory file is not created.

Note

Start with a crude calculation to get a first overview of the dynamics of a system and its possible stable or meta-stable states. Later refine parameters to ensure the quality of structural data. Such a run could look like this: Use a high starting value for Temperature initial, a low value for Temperature end, low precision and soft potentials (_s), a single k-point (gamma), Real space integration and limit SCF-convergence to 1.0E-3 or 1.0E-4 eV.

Canonical (nVT): Molecular dynamics at a constant number of particles, n, constant volume, V, and constant temperature, T, making use of the Nosé thermostat. A description of the theoretical background is given by Nosé [7] and references therein. In addition to the other parameters mentioned above the system dependent Nosé mass should be chosen.

Nose mass: The Nosé mass controls the frequency of temperature oscillations during the simulation and should be chosen such that the temperature fluctuation occurs at about the same frequency as typical phonon modes of the system. If the Nosé mass is not set explicitly, it will be chosen such that temperature fluctuates with a period of 40 time steps. The approximate frequency of the temperature fluctuations induced by the thermostat is reported in OUTCAR.out.

Isothermal-isobaric (nPT): Parinello-Rahman molecular dynamics at a constant number of particles, n, constant pressure, P, and constant temperature, T, making use of the Langevin thermostat [8], [9] . The thermostat requires a friction coefficient to be set for each atomic species as well as for the lattice degrees of freedom. These coefficients are automatically set to 50 ps^-1 for the atoms and 10 ps^-1 for the lattice, and can be modified by specifying the LANGEVIN_GAMMA and LANGEVIN_GAMMA_L tags for atom and lattice related parameters, respectively, in the Add to Input Tab.

Parinello-Rahman mass: The Parinello-Rahman mass is a fictitious mass parameter for lattice degrees of freedom (Parinello-Rahman barostat [8], [9] ) in units of a mass (AMU). Chosen too large results in a very slow variation of lattice degrees of freedom and inefficient sampling, while chosen too small leads to large geometric changes and numerical problems. If the Parinello-Rahman mass is not specified, it is set to 1000 AMU.

Continuation of job: Previous molecular dynamics calculations can be continued by specifying the previous job number in the entry field or by clicking the … button and selecting the job number from a browser. To continue a previous molecular dynamics simulation, the VASP user interface must be launched for the final configuration (finalconfig.sci) of the previous job to be continued. The molecular dynamics calculation then is continued properly by making use of the final configuration and the final velocities of the previous job. Note that the continuation of molecular dynamics is only possible using the same version of the VASP GUI, e.g. exclusively VASP 5.4 can be used to continue a previous job run by VASP 5.4. Compatibility of VASP versions and consistency of the active system and the final configuration system is checked and warnings are issued in case of discrepancies.

2.10.5.1.4. Time-dependent hybrid / DFT:

Time-dependent Hartree-Fock, hybrid functional or DFT calculations solving the Casida equations [10], as well as solving the Bethe-Salpeter equation on top of GW quasiparticle calculations (see next section) are two different approaches to obtain very accurately the response functions and thereby the optical spectra including excitonic effects. Electron-hole interactions can have a dramatic effect on the optical spectra, adding additional excitonic absorption bands.

These calculations are quite CPU time and memory intensive, and the provided parameters can be used to fine-tune computational efficiency and accuracy. Some testing might be necessary. If entry fields are left empty the VASP defaults will be applied. All parameters are discussed in the section on the GW calculations below.

2.10.5.1.5. Quasiparticle Spectra (GW)

Quasiparticles are excited states of a many-body system with energies relatively close to the ground state. The so-called GW approach truncates a series expansion of the electronic self-energy in terms of the single-particle Greens function G and the screened Coulomb interaction W, and is known as one of the most accurate approaches to calculate the excited states of a solid-state system [11], [12]. The method requires knowledge of the fully frequency-dependent dielectric function, which is obtained from the electronic wave functions from a density functional or hybrid functional-based electronic structure calculation.

Selecting Quasiparticle Spectra (GW) from the Type of calculation pull-down menu, it is necessary to define a standard DFT or hybrid functional simulation including relevant parameters in other panels. However, self-energies can only be calculated for the k-mesh as specified in the SCF panel, and properties such as the band structure, density of states, optical spectra, and response tensors requiring a different set of k-points are not available. Consequently, it is not possible to obtain a quasi-particle band structure for any chosen path through the Brillouin zone by the current version of VASP.

Three different types of GW calculation are available:

quasiparticle shifts: This evaluates the quasiparticle shifts which need to be applied to obtain correct excitation energies from electronic states derived from DFT or hybrid-functionals.

excitonic effects (Bethe-Salpeter): This adds on top of the quasiparticle spectra obtained by the above calculation type a further step solving the Bethe-Salpeter equation for electron-hole pairs to obtain the response functions and optical spectra including excitonic effects [13].

frequency dependent self-energy: This evaluates the frequency-dependent self-energies (eigenvalues of the self-energy operator) bypassing the calculation of quasiparticle shifts by the two calculation types above.

Several different procedures established for the GW method are selected by the Update pull-down menu. This specifies the extent of self-consistency applied in the procedure, i.e. whether or not eigenvalues and/or wave functions are updated in the SCF procedures for G and/or W. The resulting techniques are rather different in terms of accuracy and robustness, and are briefly discussed below:

response functions (optical spectra): This calculates the frequency dependent response functions, i.e. the dielectric function, bypassing any additional efforts to obtain quasiparticle shifts. It is emphasized, that the imaginary and real part of the frequency-dependent dielectric functions are always determined in the course of a GW calculation. It can be assembled by looking for the string “dielectric constant” (two blanks between the words) in GW_OUTCAR.out. Two different data sets can be found: the first one is the head of the microscopic dielectric matrix not including local field effects, the second data set is the inverse dielectric matrix including local field effects.

Local field effects can be treated on the Hartree level based on the random phase approximation (RPA) or including changes of the exchange-correlation potential on the density functional level (DFT). The DFT level is the default, which can be changed to RPA by adding the line “LRPA = .TRUE.” into the Add to Input panel. These options are available also for DFT and hybrid functional-based optical properties calculations.

no updates, perturbative (G0W0): This performs a single perturbative GW step without updating the eigenvalues and wave functions of the preceding DFT or hybrid functional calculation [14].

Update of eigenvalues for G (GW0): This applies a partially self-consistent GW algorithm with an iterative update of eigenvalues for the calculation of the Greens function only. The wave functions of the preceding ground state calculation based on DFT or hybrid functionals are not updated. As a default, 4 update steps are applied. This is the default update option because for most cases this procedure yields results closest to experimental data with much less computational demands than the below options [15].

Update of eigenvalues for G and W (GW): A partially self-consistent GW algorithm iteratively updating the eigenvalues for the calculation of the Greens functions G and the screened Coulomb potential W is applied. In update options, the wave functions of the preceding ground-state calculations are used without update. This option tends to decrease the agreement to experimental data on band gaps in most cases.

Update of wave functions & eigenvalues for G: A self-consistent GW algorithm is applied with a full update of wave functions and eigenvalues for the calculation of the Greens function. Although the computational efforts are drastically increased the results are typically less satisfactory as those obtained from the GW0 approach. Furthermore, update of the wave functions in general decreases the level of the robustness of the algorithm.

Update of wave functions & eigenvalues for G and W: A fully self-consistent GW algorithm updating wave functions and eigenvalues for the calculation of both G and W improves the agreement to experimental data close to or above the level achieved by the GW0 method, however, with much larger computational efforts and the considerable probability of computational issues. The main advantage can be gained by including so-called vertex corrections, taking into account electron-hole interactions [16]. Vertex corrections, however, are still an undocumented feature of the VASP code and are very difficult and tedious to perform.

GW calculations are extremely demanding in terms of CPU time and memory allocation, and it might frequently be necessary to fine-tune specific parameters affecting computational efficiency to make these calculations feasible. Of course, one has to trade efficiency against accuracy and some testing will be indispensable. The following parameters allow for the fine-tuning of the GW process. If entry fields are left empty the VASP defaults will be applied, which may not be suitable for the system to be studied.

Restore charge density: This specifies how the all-electron charge density is restored on the plane wave grid and influences the accuracy of the eigenvalues. For hybrid functional, screened exchange and Hartree-Fock calculations the moments of the charge density only are restored as a default. This option is made available for GW calculations by selecting moments only. For GW calculations, in particular for systems with localized electronic states, it is recommended to restore the full shape of the all-electron density on the plane wave grid up to a certain angular momentum lmax. For first and second row elements the all-electron charge density with full shape up to lmax = 2 should be restored, for transition metal elements the option full shape up to lmax = 4, and for f-electron systems the option full shape up to lmax = 6 is more appropriate.

Number of update iterations: This specifies the number of iterations performed for achieving self -consistency for any of those Update options above requiring updates of eigenvalues and/or wave functions for G and/or W. The default value of 4 update iterations has proved sufficient for many cases.

Cutoff for response functions: This defines the basis set for the response functions in the same way as the plane wave cutoff defines the basis set for the wave functions (see Plane wave cutoff in section c, General Setup). If this parameter is not set explicitly in this entry field, it is set to 2/3 of the plane wave cutoff, which yields reasonable accuracy at a moderate computational cost. Values between 150 and 200 eV are found to be sufficient, for some cases even 100 eV may be suitable. Specification of a lower cutoff speeds up calculations and reduces memory demands substantially. For convergence tests it is recommended to increase the Plane wave cutoff keeping Cutoff for response functions at a constant ratio.

Number of bands: The number of bands used for the calculation of response functions and optical spectra in GW must include a large number of unoccupied bands. If not set explicitly, the default number of bands as shown above, the entry field is applied, which is a much larger number then the default number of bands used for standard SCF calculations. The specified number may become slightly increased for parallel runs in order to obtain a multiple of the number of processors operating in parallel.

Number of quasiparticle bands: This is the number of bands for which quasiparticle shifts should be calculated by GW. As a default the value for standard SCF calculations as shown above the entry field is used, which may involve many more bands than are of interest, thereby increasing drastically the computational efforts. Always choose this value such that only the states of interest are covered.

Number of frequency points: This specifies the number of frequency grid points for the evaluation of frequency-dependent functions in GW. This number should be chosen around 50-100. For parallel runs, the number should be a multiple of the number of compute nodes for maximum efficiency. For quick and less memory demanding calculations values around 20-30 are sufficient. With this setting, however, errors of the order of 20-50 meV for the gap and 100-200 meV for the bottom of the conduction band must be expected. It is not recommended to increase this value beyond 100 for a k-point sampling of 4x4x4 k-points/atom because the joint density of states and the self-energy tend to exhibit a spurious fine structure related to the finite k-point grid. This fine structure is smoothed by a smaller number of frequency points or by more k-points. For a 6x6x6 k-points/atom grid the number of frequency points can usually be increased to 200-300 without noticing problems associated with this kind of noise. It is noted that this parameter does not influence critically the CPU time demands.

Complex shift parameter: The small complex shift should be chosen depending on the number of frequency points, i.e. for less dense frequency grids the shift parameter should be increased accordingly. As default behavior the complex shift parameter is not set in the entry field and VASP determines the value such that the calculations are converged to about 10 meV with respect to the number of frequency points. This means that for constant complex shift parameter, the quasiparticle shifts should not change by more than 10 meV, if the number of frequency points is increased. The parameter should be at least as large as the grid spacing at low frequencies, if chosen smaller the quasiparticle energies might show erratic behavior.

For the excitonic effects (Bethe-Salpeter) type of GW calculation two additional parameter occur at the bottom of the subpanel which can be used to tune computational demands and accuracy:

Number of occupied bands: determines how many occupied orbitals, counted from the Fermi energy, are included in the Bethe-Salpeter calculations

Number of unoccupied bands: determines how many unoccupied (virtual) orbitals just above the Fermi level are included in the Bethe-Salpeter calculations

Note that compute time for Bethe-Salpeter calculations increases with the third power of the number of occupied/unoccupied bands, whereas the memory demands increase quadratically. Default values are shown above the entry field. For highly accurate results, the Number of unoccupied bands often needs to be increased, whereas for large systems one is often forced to reduce both values too much smaller numbers.

2.10.5.1.6. Accurate Energy (ACFDT-RPA)

This calculation type evaluates a very accurate total energy for the input structure (without geometry optimization) as a sum of the exact exchange energy and the correlation energy within the random phase approximation (RPA) utilizing the adiabatic connection fluctuation dissipation theorem (ACFDT) [17], [18].

These calculations are CPU time and memory intensive, and the provided parameters can be used to fine-tune computational efficiency and accuracy. Some testing might be necessary. If entry fields are left empty the VASP defaults will be applied. Most parameters are discussed in the section on GW calculations above, but there are a few special aspects:

Number of frequency points: much smaller default value of 12 is needed as compared to GW and time-dependent DFT/HF calculations aiming at the response functions on that grid. For large gap systems one might obtain good convergence alcreated using 8 points, whereas for metals up to 24 frequency points are sometimes necessary, in particular, for large unit cells.

Run the ACFDT-RPA algorithm for metallic systems: this makes sure that exchange and correlation energy are evaluated on the same k-point grid, long-wavelength contributions from the polarizability are not considered, and a correction energy for the exact exchange energy related to partial occupancies is added.

2.10.5.1.7. MT-Elastic Properties

The optional MT- Elastic Properties module is tightly integrated into the VASP interface.

The actual determination occurs in three steps:

• MT analyzes the cell symmetry and determines the required directions of strain to derive elastic constants.
• The JobServer performs a VASP calculation for each of the strained cells.
• MT analyzes the results and computes the elastic constants (elastic constants and compliances matrices together with their Eigenvalues and Eigenvectors). In addition, elastic moduli and sound velocity, as well as derived thermodynamic properties (within the framework of the Debye model) are reported.

Note

For accurate results, optimize the initial cell first to minimize residual forces and stresses before running MT-Elastic Properties. Note that elastic constants depend considerably on the lattice parameters (volume) of the system. You should enforce a much tighter Convergence than the default value of 0.02 eV/\({\mathring{\mathrm{A}}}\).

The most important parameter is the amount of strain: The resulting forces should be as large as possible (to get a better signal to noise ratio) without leaving the confines of the elastic regime.

Strains: a list of strains in fractions of the unit cell. The value 0.005 refers to 0.5% strain. You can add more than one strain amount, separated by spaces, to get better results, when the material is anisotropic, e.g. 0.005 0.01.

Start from wave functions of unstrained structures: In general, this option is not recommended, since VASP may have difficulties to adapt wave functions to distorted cells resulting in numerical instability. This option may be useful to retain magnetic ordering.

Relax atom positions of strained structures: Straining the structure can break the symmetry of the cell and leaves one or more degrees of freedom for atoms to relax. Check this option if you want to optimize the atom positions of these cells. For comparison to experimental data, it is mandatory to optimize all atomic degrees of freedom for all distorted unit cells.

If Relax atom positions of strained structures is checked, entry fields for two parameters guiding the optimization of atom positions become available, which coincide with the same parameters of the Optimization subpanel (see section ii, Structure Optimization):

Convergence: A much tighter convergence criterion than the default value of 0.02 eV/\({\mathring{\mathrm{A}}}\) is recommended. For instance, 0.002 eV/\({\mathring{\mathrm{A}}}\) would be a suitable choice.

Maximum number of steps: Sets the maximum number of geometry steps.

Note

To obtain accurate elastic constants, it is recommended to use fine k-point sampling in reciprocal space (e.g. the spacing of k-points of 0.2/\({\mathring{\mathrm{A}}}\) or below for metallic systems) in the k-mesh section of the SCF tab. Furthermore, the linear tetrahedron method was found to provide reliable results and by checking Extrafine augmentation grid for accurate forces in the Advanced/Restart panel the accuracy of the optimizations maybe enhanced.

2.10.5.2. Properties

Check items in the Properties frame to calculate and write out the related properties following a VASP calculation as selected by the Type of calculation pull-down menu (see previous section Type of calculation).

Note

The above options are available for all types of calculations. However, depending on your target it might be more efficient to run two independent steps to get specific properties.

Example: You are running a structure optimization and you would like to determine the total energy, the band structure, and the DOS for the resulting system. During your structure optimization, the shape and volume of your starting unit cell may change significantly (>3%). In this case, results for total energies and related properties will improve if you run the optimization first, reload the optimized system and then start a new VASP job to get the desired properties. The reason is that some VASP parameters implicitly depend on the input geometry, i.e. the cell shape and size. For consistency, MedeA does not change these parameters in the course of a set of tasks within a given job.

(Pseudo, difference, spin) charge density: The electronic pseudo charge density and several derived data are provided. This includes:

• The raw data for the pseudo charge density as provided by VASP (data written to CHGCAR). This file can be used for the restart of related calculations, e.g. to save CPU time or stabilize magnetic states. Restart is handled by the options of the Initial conditions and restart frame of the Advanced/Restart panel. This charge density is called pseudo charge density because of unphysical shapes within spheres of the core radii (so-called depletion radii).
• Previous to MedeA version 2.12, i.e. before total charge densities became available (see below), the so-called valence charge density was provided, with unphysical parts within the core radii replaced by a steep and radially symmetric charge density rise approaching the positions of the nuclei (data written to ValenceChargeDensity.data). As of MedeA version 2.12 this data file is not supported anymore; however, previously generated ones can still be viewed and analyzed.
• The pseudo charge density obtained from the superposition of atomic charge densities placed at their lattice positions (data written to ATOMS..CHGCAR).
• The difference charge density between the self-consistent pseudo charge density (CHGCAR) and the superposition of atomic charge densities (ATOMS..CHGCAR), which sometimes is referred to as deformation charge density (data written to DifferenceChargeDensity.data). Unphysical parts within spheres of core radii cancel out.
• The magnetization density, i.e. the difference between spin-up and spin-down pseudo charge densities, is provided for all spin-polarized calculations (require data written to CHGCAR). Unphysical parts within spheres of core radii cancel out.
• The magnetization density in Cartesian x, y and z direction is provided for all non-collinear and spin-orbit magnetic calculations.

Total local potential: the Coulomb potential, excluding the exchange-correlation potential (data written to LOCPOT)

Electron localization function: The ELF is a particular way to analyze the wave functions to understand chemical bonding (data written to ELFCAR)

Wave functions: The electronic wave functions. Saving the wave functions requires a lot of hard-disk space, but is very useful when planning to later restart a calculation (data written to WAVECAR.txt, which is a translation of the WAVECAR file directly written by VASP into a machine independent ASCI file).

Electric field gradients: The electric field gradients at the positions of the nuclei as well as the quadrupolar coupling constants as measured in nuclear magnetic resonance experiments. The settings for the nuclear quadrupole moments of the nuclei [19] (the QUAD_EFG tag in INCAR) can be viewed from the Preview Input Tab, which also provides alternative settings if more than one isotope with nuclear quadrupole moments are available:

LEFG = .TRUE.
QUAD_EFG = 0 146.6 -25.58 33.27 2.860 # Al-27 O-17 C-11 H-2
# Nuclear electric quadrupole moments of other isotopes are not
available in the database.

Alternative settings can be applied by adding the appropriate QUAD_EFG tag to the INCAR file from the Add to Input Tab of the VASP GUI. All results are collected in Job.out.

Hyperfine parameters: The hyperfine parameters describing the interaction between the spin of the nuclei and the electronic spin density. This property requires a magnetic setup (see section 4 on definitions for the interaction). If a non-magnetic setup is specified the checkbox is inactive and its setting is ignored.

Work function (surfaces only): Use for surface models only. The energy required moving an electron from the top of the valence band to infinity. The work functions are reported in Job.out

So far, any of the above properties are provided from the simulation as specified by the Type of calculation. Except for the evaluation of the superposed atomic charge density, this does not require extra tasks. The following properties, however, require additional tasks, which may considerably increase computational demands.

(Total, valence) charge, Bader analysis: The total electronic charge density and several derived data are provided. This includes:

• The total charge density including valence and core electrons, thus being a true physical observable (data written to CHARGES CHGCAR_total). This is the sum of the total valence charge density and the total core charge density.
• The total valence charge density (data written to CHARGES_AECCAR2).
• The core charge density (data written to CHARGES_AECCAR0).
• The total valence charge density obtained from the superposition of total atomic charge densities placed at their lattice positions (data written to CHARGES_AECCAR1).
• The difference charge density between the self-consistent total valence charge density (CHARGES_AECCAR2) and the superposition of total atomic charge densities (CHARGES_AECCAR1), which sometimes is referred to as deformation charge density (data written to CHARGES CHGCAR_totaldiff).
• The Bader charge analysis [20] is summarized in Job.out, providing charges, charge transfer with respect to atoms, and the result of the Bader volume decomposition.

Band structure: Dispersion relation E(k), i.e. the electronic energy as a function of momentum along a path through the Brillouin zone (data written to BandStructure.data). Further parameters can be set in the Band Structure panel.

Electronic Density of states (DOS): The distribution of the number of electronic states as a function of the energy (data written to DOSCAR, DensityOfStates.data). Further parameters can be set in the DOS/Optic/Tensors panel.

Optical spectra: The frequency dependent optical spectra, such as real and imaginary part of the dielectric function and conductivity, furthermore reflectivity, adsorption, and refractory index (data written to OpticalSpectra.data). If this property is requested, the DOS will be provided as well. Further parameters can be set in the DOS/Optic/Tensors panel. In particular, it is strongly recommended to avoid the tetrahedron method, since it may cause errors. Gaussian smearing or one of the integration techniques other than the tetrahedron method are recommended for optical spectra. The density of states could be run separately with the tetrahedron method if needed.

Zone center phonons: A finite differences approach is applied to evaluate phonon frequencies at the \(\Gamma\) point. Opposite to the linear response approach, this does work for hybrid functional as well.

Response tensors: This runs linear response (density functional perturbation theory) to obtain the dielectric tensor (low and high frequency limit), the piezoelectric tensor, Born effective charges, and zone center phonon frequencies. The tensors are reported in Job.out and refer to the unit cell and atoms as defined internally in VASP (POSCAR). Further parameters can be set in the DOS/Optic/Tensors panel. Linear response is only available for density functionals, but not for methods involving non-local Hartree-Fock exchange (the checkbox is grayed out if such a functional is chosen).

NMR: chemical shifts: This runs linear response calculations to obtain chemical shifts [21], [22] as measured in nuclear magnetic resonance (NMR) experiments. The results are listed in Job.out. Note that only relative chemical shift data of atoms of the same element on different sites are physically meaningful.

The following graphical visualization options for these properties are available in MedeA (see also section II. H. Analysis of Results):

• (VASP) Trajectories, Trajectories, Gibbs Trajectories
• Band Structure
• Density of States
• Optical Spectra
• Total Charge Density
• Total Valence Charge Density
• Difference Charge Density
• Magnetization Density
• Total Charge Density and Total Valence Charge Density
• Electron Localization Function
• Total Local Potential

Suitable for molecules in an otherwise empty box or for surface models, an implicit solvation model can be applied (VASPsol [23] ) simulating the effect of a solvent on the geometry, total energy, and dynamics of the system.

Apply solvation model: Switches on the implicit solvation model and providing access to

Solvent dielectric constant: The solvent’s bulk dielectric constant, which may be known from experimental data. The default value of 78.4 applies to water and covers mostly the electrostatic interactions. Other default parameters defining the dielectric cavity, and capturing non-electrostatic effects of cavitation, dispersion, and repulsion are optimized for water and can be modified for other solvents from the Add to Input Tab, as explained in detail by the context-sensitive Help text (however, optimized values are currently only available for water).

There is an external condition that can be set for all simulations:

External pressure: The hydrostatic pressure can be specified in units of GPa, a positive value indicating compression. To have the system responding to the external pressure, two subsequent geometry optimizations of all structural degrees of freedom (cell volume, cell shape and atom positions) are recommended, before further simulations are performed.

For the calculation of charged complexes or defects, the charge state of the system can be set for all simulations:

Charge state: The charge of the system can be specified in units of the elementary positive charge e. A negative charge state means excess electronic charges, i.e. extra electrons added to the system, whereas a positive charge state means a lack of electrons. To conserve overall charge neutrality, a compensating homogeneous background charge is assumed.

Besides, an external electric field can be applied:

External electrostatic field: This choice allows whether or not to apply an external electrostatic field in a given Cartesian direction. This is a valid option only for slab and molecular systems. Restart from wave functions generated for the system without an electric field is recommended.

Electric field strength: In case an external electrostatic should be applied this entry enables to specify the field strength in units of eV/Ang in the chosen Cartesian direction.

2.10.5.3. Interaction

The parameters of this panel define the level of theory applied for the quantum-physical description of electron-nuclei and electron-electron interactions. The Functional menu in combination with DFT exchange-correlation allows you to select the degree of locality, spanning the range between the local density approximation of density functional theory (DFT) up to pure Hartree-Fock with several intermediate (hybrid) steps. Van der Waals applies a force field-based correction. Magnetism provides several options for including magnetic effects, ranging from non-polarized, and polarized up to non-collinear spin-orbit coupled Hamiltonians.

Parameters in this panel affect all calculation types and property calculations. In particular the choices for the Functional and the DFT exchange-correlation and the detailed choice of potentials for each element of the system provided in the Functional/Potential panel critically determine the absolute values of calculated total energies. This means that total energies can only be compared if these parameters are chosen consistently. To a large extent this is also true for parameters such as Magnetism, Planewave cutoff, and to a smaller extent for the Precision, as provided in the General Setup frame below. On the other hand, measurable properties such as optimized geometries, elastic constants, heats of formation and thermodynamic properties rely in most cases much less critically on consistent and transparent use of these settings. It should be emphasized, however, that it is a good and strongly recommended practice to keep these parameters as consistent as possible for a given study.

The specific options of the Interaction frame are:

Functional: This provides a fundamental choice of how to treat on a quantum-physical basis the interactions between the electrons and between the electrons and the nuclei. The options are:

Density functional - The interactions are determined from the total density. The exchange correlation may be calculated from the local density only (local density approximations, LDA) or, also, from the gradient of the local density and the knowledge about the rules guiding the shape of the exchange-correlation hole (generalized gradient approximation, GGA). These options are provided by the

DFT exchange-correlation choice, as discussed below.

Van der Waals density functional - The van der Waals density functionals (optB86b-vdW, optB88-vdW, optPBE-vdW, BEEF-vdW, rev-vdW-DF2, rPW86-vdW2, revPBE-vdW, SCAN + rVV10) approximately take into account the dispersive forces and van der Waals interactions. These approaches do not rely on empirical forcefields (which can be selected from the option Van der Waals below) and are true first principles techniques. The exchange correlation is predefined with the chosen functional, therefore the option DFT exchange-correlation is inactive and its setting is ignored. Furthermore, also the Van der Waals forcefields are not applicable and the option Van der Waals is inactive, therefore. If this option is selected, an additional menu Type of van der Waals functional appears, which allows you to specify the specific functional.

Meta-GGA - The meta-GGAs (revTPSS, TPSS, SCAN, MS2, MS1, MS0, M06-L, MBJLDA) in addition to the density make use of the kinetic energy density for more accurate energies and structures. The exchange correlation is predefined with the chosen functional, therefore the option DFT exchange-correlation is inactive and its setting is ignored. If this option is selected, an additional menu Type of meta-GGA appears, which allows you to specify the specific functional.

Hybrid functional - The interactions are evaluated from a mix of the density functional local or semi-local exchange and correlation and the exact non-local Hartree-Fock exchange. If this option is selected, an additional menu Type of hybrid functional appears, which allows you to specify different realizations of this concept.

Screened exchange - The correlation is treated as density functional and the screened exchange contribution - as described as density functional - is replaced by the Thomas-Fermi screened non-local Hartree-Fock exchange.

Hartree-Fock - The interaction is treated via exact non-local Hartree-Fock exchange only without screening and correlation is not applied.

DFT exchange-correlation: Different approximations to the exchange-correlation part of the density functional can be made. The available options are:

LDA - The local density approximation as parameterized by Perdew & Zunger [24] with electron correlation obtained from Quantum Monte Carlo simulations of Ceperley & Alder [25]

GGA-AM05 - The generalized gradient approximation (GGA) after Armiento & Mattsson [26]

GGA-PBEsol - The GGA tuned in particular for properties of solids [27]

GGA-PBE - The general purpose standard GGA after Perdew, Burke & Ernzerhof [28]

GGA-rPBE - The GGA tuned for adsorption and surface properties [29] by Hammer, Hansen & Norskov

GGA-BLYP - The GGA after Becke [30]

Type of van der Waals functional: Different types of van der Waals functionals have been implemented in VASP, and can be accessed by this pull-down menu:

optB86b-vdW [31] - The optimized van der Waals functional based on the Becke 86 exchange functional, which tend to exhibit the smallest errors for most systems investigated.

optB88-vdW - The optimized van der Waals functional based on the Becke 88 exchange functional with accuracies comparable to optB86b-vdW for most systems investigated.

optPBE-vdW - The van der Waals functional based on the PBE exchange functional with optimized enhancement factors.

BEEF-vdW [32] - The Bayesian error estimation functional developed by J. Wellendorff and co-workers.

rev-vdW-DF2 [33] - The van der Waals functional developed by I. Hamada.

rPW86-vdW2 [34] - The van der Waals functional based on the Perdew-Wang 86 exchange functional.

revPBE-vdW [35] - The original van der Waals density functional of Dion et al. making use of the rPBE exchange approximation, which tends to larger errors than the optimized functionals above.

SCAN + rVV10 [36] - is supplementing the strongly constrained and appropriately normed (SCAN) meta-generalized gradient approximation for short- and intermediate-range interactions with the long-range vdW interaction from rVV10, the revised Vydrov-van Voorhis nonlocal correlation functional [37], [38].

For comparison of the performance of the various types of van der Waals density functionals, see Klimes et al..

Type of meta-GGA: Different types of meta-GGA functionals have been implemented in VASP, and can be accessed by this pull-down menu:

revTPSS - The revised Tao-Perdew-Staroverov-Scuseria [39] , [40] functional improves surface energies and atomization energies as well as lattice parameters combining the advantages of the TPSS meta-GGA and the PBEsol density functional.

TPSS - The original Tao-Perdew-Staroverov-Scuseria [41] functional provides improved surface energies and atomization energies but only minor improvement for lattice parameters. This is achieved by respecting two paradigms, the uniform electron gas in condensed matter physics and the hydrogen atom for chemistry.

SCAN [42] - The Strongly Constrained and Appropriately Normed meta-GGA functional fulfills all known constraints that the exact density functional must fulfill. There are indications that this functional is superior to most gradient corrected functionals.

MS2, MS1, MS0 [43] , [44] - The series of so-called Made Simple functionals are believed to improve the description of noncovalent interactions over PBE, TPSS and revTPSS, but not over M06-L.

M06-L [45] - This meta-GGA functional is constructed to satisfy the uniform electron gas limit and is fitted to molecular data to perform well for main-group and transition metal chemistry.

MBJLDA [46] - The Becke-Johnson exchange potential combined with the LSDA correlation yields band gaps with an accuracy comparable to hybrid functional or GW calculations, however, with computational demands comparable to standard density functional computations.

Attention

Since MBJLDA applies the LSDA exchange-correlation energy instead of an appropriate one, it is a potential-only functional. As a consequence MBJLDA is not self-consistent with respect to the total energy, Hellman-Feynman forces cannot be computed and geometries cannot be optimized. The functional is aimed at the calculation of band structures, in particular band gaps, densities of states and optical properties and geometries should be optimized by a different functional in a separate job.

Note

MBJLDA calculations for surfaces tend to diverge, because of the functional becomes unstable in vacuum.

Type of hybrid functional: Different types of hybrid functionals have been implemented in VASP, and can be accessed by this pull-down menu:

HSE06 - The hybrid functional developed by Heyd, Scuseria & Ernzerhof [47] exhibits good convergence with respect to k-point sampling in the Brillouin zone and therefore is the recommended default approach.

PBE0 - The standard hybrid functional for solids developed by Ernzerhof & Scuseria [48] and Adamo & Barone [49].

B3LYP - The standard hybrid functional for molecules after Becke [50], Lee, Yang & Parr [51], Vosko, Wilk & Nusair [52] , Stephens, Devlin, Chabalowski & Frisch [53].

Note

The hybrid functionals HSE06 and PBE0 chosen as Type of hybrid functional has been developed with a density functional approximation for exchange-correlation as described by GGA-PBE [47]. There is, however, no reason not to use other approximations as offered by the DFT exchange-correlation menu. Even LDA is a sensitive choice to be combined with a hybrid functional. The MedeA implementation, therefore, allows any combination of HSE06 and PBE0 with all of the approximations for exchange-correlation. On the other hand, the B3LYP functional works only together with the GGA-BLYP approximation.

It is noted, that calculations based on non-local Hartree-Fock exchange become more efficient or may even be feasible only if starting from converged wave functions obtained within a density-functional. Therefore, protocols are implemented in MedeA enabling full automation of such multi-step procedures. This is discussed in more detail in 2. Functionals).

Finally, it is noted that the k-point sampling requires special attention when methods based on non-local Hartree-Fock exchange are applied. Default settings can be reviewed in the SCF panel for the calculation types (see section 6. The SCF panel), and in the DOS/Optic/Tensors panel for the density of states, optical properties and response tensors (see section II.G. 7. The DOS/Optic/Tensors panel). The contents of these panels are depending on whether density functional methods or approaches based on Hartree-Fock exchange are applied.

Van der Waals: This adds a Van der Waals contribution to the interactions specified above. Van der Waals contributions are provided only for selected functionals. The choices are:

None - No van der Waals contribution is added.

DFT-D3 zero-damping [54] - The interactions are added by a pairwise additive forcefield with atom-pairwise specific dispersion coefficients and cutoff radii that are both computed from first principles, as well as eighth-order dispersion terms. System and geometry specific information is used by fractional coordination numbers. A standard zero-damping formula is applied. Forcefield parameters and scaling factors are available for PBEsol, PBE, rPBE and BLYP, the meta-GGAs TPSS and M06-L, and hybrid functionals PBE0 and B3LYP, covering all elements up to Pu.

DFT-D3 BJ-damping [55] - The same scheme for computation of the dispersion coefficients are applied than above. However, an alternative rational damping to finite values for small interatomic distances according to Becke and Johnson (BJ-damping) is used, which avoids repulsive interatomic forces at shorter distances. Forcefield parameters and scaling factors are available for PBEsol, PBE, rPBE and BLYP, the meta-GGA TPSS, and hybrid functionals HSE06, PBE0 and B3LYP, covering all elements up to Pu.

Tkatchenko-Scheffler [56] - The interactions are added by pairwise additive semi-empirical forcefield with C6 parameters calculated from the self-consistent DFT charge density. As of MedeA 2.18 for improved numerical stability the standard approach is applied rather than the Hirshfeld partitioning. These forcefield parameters and scaling factors are available for PBE , the meta-GGAs TPSS and M06-L, and hybrid functionals HSE06, PBE0 and B3LYP.

Tkatchenko-Scheffler + SCS [57] - The Tkatchenko-Scheffler van der Waals method is combined with the self-consistent screening equation of classical electrodynamics to describe microscopically the frequency-dependent polarizability of finite-gap molecules and solids (Phys. Rev. Lett. 108, 236402 (2012); Phys. Rev. B. 87, 064110 (2013)). Forcefield parameters and scaling factors are available for PBE.

Many-body dispersion energy [58], [59] - is based on the random phase expression for the correlation energy, whereby the response function is approximated by a sum of atomic contributions represented by quantum harmonic oscillators.

DFT-dDsC dispersion correction [60] - The approach is closely related to DFT-D2, but is charge density dependent. A sufficiently dense FFT grid, and Accurate Precision must be therefore recommended for such calculations.

DFT-D2 forcefield (Grimme) [61] - The interactions are added by pairwise additive semi-empirical forcefield as suggested by Stefan Grimme. Forcefield parameters and scaling factors are available for PBE, BLYP and B3LYP functionals.

Magnetism: Magnetism plays an important role in the structure of some metals and many molecules, oxides, and oxide surfaces. Options for capturing magnetic effects are

Defined by model - MedeA allows you to set magnetic moments for each atom in the structure window via the context menu obtained by a right-mouse-click on an atom (Atom/Magnetic Moments…), or in the Molecular Spreadsheet. (Note, that also the atomic mass of each atom can be modified through this panel or the Molecular Spreadsheet, e.g. for studying isotope effects in dynamics or vibrational analysis.) If magnetic moments are set for one or more atoms of the system, a spin-polarized VASP calculation is performed and these moments will be used to set up the initial magnetic structure. If no magnetic moments are set on any of the atoms, a non-magnetic VASP calculation is performed.

Non-magnetic - No magnetism will be considered

Spin-polarized - A magnetic (spin-polarized) calculation is performed allowing for scalar spin-up and spin-down magnetic moments for each atom (ferromagnetic, antiferromagnetic, or antiferromagnetic spin configurations may be reached after self-consistence).

Non-collinear magnetic - A magnetic calculation is performed that accounts for non-collinear spin configurations, i.e. atoms are allowed to develop spin components in any direction of space during self-consistence cycling.

Spin-orbit magnetic - A fully relativistic magnetic calculation including spin-orbit coupling is performed (Dirac equation). The spin-quantization axis can be specified in the Advanced/Restart panel, if a spin-orbit magnetic calculation is specified (otherwise these entry fields are hidden). Initial magnetic moments are set for atoms in the model are assumed to be oriented in the Cartesian z-direction.

Spin-orbit coupling and non-collinear magnetic calculations can be fairly time-consuming compared to spin-polarized and non-magnetic calculations.

2.10.5.4. General Setup

Parameters in this panel affect all calculation types and property calculations. In particular, the potentials as defined by the Potential menu and the detailed choice for each element of the system provided in the Potentials panel, critically determine the absolute values of calculated total energies. This means that total energies can only be compared if potentials are chosen consistently. To a large extent, this is also true for parameters such as Magnetism and Planewave cutoff, and to a smaller extent for the Precision. On the other hand, measurable properties such as optimized geometries, elastic constants, heats of formation and thermodynamic properties rely in most cases much less critically on consistent and transparent use of these settings. It should be emphasized, however, that it is a good and strongly recommended practice to keep these parameters as consistent as possible for a given study.

The default settings for these parameters are carefully chosen and tested to provide reliable results for most systems without being computationally too demanding.

Precision: Influences several internal parameters such as plane-wave cutoff (basis set size), FFT (Fast Fourier Transform) mesh (normal mesh and fine mesh), fine grid Fourier mesh, and integration mesh for real space projectors

Normal - For standard calculations

Low - For crude molecular dynamics and “first guesses”

Accurate - For precise energies and forces, where the lattice parameters remain unchanged

Standard 500 - Overall high precision cutoff, used for MedeA’s reference heat of formation energies. To entirely comply with the Standard 500 settings for accurate formation energies, the density of the k-mesh needs to be increased by selecting a Spacing of k-points of 0.2 1/\({\mathring{\mathrm{A}}}\) in the SCF panel, in addition. This setting is applicable to all types of calculations including cell parameter optimization, but maybe unnecessarily demanding in terms of CPU resources.

Increase plane wave cutoff (cell optimizations): This is required for cell volume and/or cell shape optimizations, if Standard 500 is computationally too demanding. This option increases the cutoff by 30 %. An increased plane wave cutoff is mandatory for cell optimizations because stress tensors are converging much slower with the number of plane waves than atomic forces.

Consider using extra fine augmentation grid (Advanced/Restart panel) for cases where extremely accurate forces are needed (difficult MT and Phonon calculations)

Plane wave cutoff: Defines the precision (=size) of the plane wave basis set. The default value determined by the current settings is shown above the entry field. You may enter a specific cutoff into the entry field, e.g. to be consistent with other calculations.

Note

It is emphasized that total energies tend to slowly converge with increasing plane wave cutoff. Therefore it is highly recommended to compare total energies of different calculations only if almost identical plane wave cutoffs have been used.

Projection: defines how the projection operators are to be evaluated

Real space (faster for larger systems, but numerically somewhat less precise)

Reciprocal space (more precise, but increasingly slower the larger the system)

For larger systems (all unit cell dim > 8 \({\mathring{\mathrm{A}}}\)) a considerable amount of time can be saved by using Projection: Real space.

Beneath the General Setup frame the pulldown menu VASP version allows to direct VASP calculations to either standard VASP executables operating with CPUs or to VASP executables suitable for GPUs. For the latter reciprocal space projection is not available, thus the Projection choice becomes grayed out.

Note

The pulldown menu VASP version is only visible if MedeA is generally enabled to run Jobs on TaskServers offering GPUs. This needs to be set by switching on Enable running Job on GPU in the Miscellaneous Tab of the Preferences window, which can be launched from the menu item File >> /Preferences….

2.10.6. The Functional/Potential Panel

The Functional/Potential panel contains two different frames, one for setting projector augmented wave potentials and one for further parameters used to tune calculations applying non-local functionals based on Hartree-Fock exchange, such as hybrid functionals and screened exchange.

2.10.6.1. Potentials

Select the Functional/Potential panel to edit the type of potentials to use in the frame to the left. In General Options you can select the approximation for exchange and correlation used for the density functional from the pull-down menu DFT exchange-correlation. This is identical to the menu appearing in the Interaction frame of the Calculation panel. The available options are the local density approximation LDA, and several flavors of the so-called semi local generalized gradient approximation, such as GGA_AM05, GGA_PBEsol, GGA_PBE, GGA_rPBE, and GGA_BLYP. Further details are discussed and references are provided in section 3. Interaction

Two different sets of the projector augmented wave (PAW) potentials are available for almost all elements of the periodic table in MedeA VASP: one set for LDA and one for GGA based calculations. In addition, for many elements several different element specific potential types are available. These can be chosen from the pull-down menus appearing for each element in the section beneath Specific Potentials per Element.

In the example shown above (GaAs), six different potentials are offered for Ga: a standard one (Ga), a hard potential (Ga_h) having a very high plane wave cut-off, and a potential incorporating the Ga semi core d-states into the set of valence states (Ga_d); In addition to these three traditional potentials, there are three further potentials with suffix GW providing a better description for unoccupied states of higher energy: a standard GW potential (Ga_GW), a GW potential incorporating the Ga semi core d-states (Ga_d_GW), and a GW potential incorporating in addition to these d-states also semi core s-states and p-states (Ga_sv_GW). Depending on your choice, the default plane wave cutoff energy is automatically updated. The terminology of different potential types is briefly summarized below:

Types of Specific Potentials per Element X:

X - standard potential

X_d - treats semi core d-states as valence states

X_pv - treats semi core p-states as valence states

X_sv - treats the semi core s-states as valence states

X_s - soft potential: soft potentials require a lower plane wave cutoff and are therefore quite fast, but less precise. They are useful for less accurate calculations, but should be avoided for short bonds and molecular systems.

X_h - hard potential: hard potentials are very precise but calculations are time-consuming due to very high plane wave cutoff energies. They are indicated for very short bonds and systems under extremely high compression.

In addition, the presence of an additional suffix GW indicates a better description of excited states, thus enhancing the accuracy for optical spectra and quasi particle calculations.

The default potentials are carefully chosen under the recommendations of the authors, providing a good compromise between accuracy and computational efficiency.

Note

In the MedeA implementation of VASP 5 the Ultrasoft Pseudo Potentials (US) are not available anymore, only Projector Augmented Wave (PAW) potentials LDA and PBE are provided. The US potentials are not supported anymore by the authors of VASP and many features of VASP 5 do not work with this class of potentials. Since the PAW recovers the correct nodal structure of electronic wave functions near the nuclei, they are more accurate than US potentials. Furthermore, the set of GGA_PW potentials based on the exchange-correlation functional of Perdew and Wang are not supported anymore by the authors of VASP and are therefore not provided in MedeA VASP 5.

2.10.6.2. Functionals

The frame to the right of the Functional/Potential panel allows you to customize and tune the functional as selected from the Functional pull-down menu, which is identical to the menu of the Calculation panel (see section 3. Interaction).

No additional settings and frames become available if Density functional is chosen.

If Hybrid functional is selected, the Type of hybrid functional can be chosen equivalently to the identical menu in the Calculation panel.

For Hartree-Fock no specific further options can be set.

If Screened exchange is chosen as Functional, two different modes for the specification of the screening term are available, i.e. it is possible to select the

Thomas Fermi screening length from:

average valence density: The screening length is automatically calculated from the average valence density as defined by the potential. As a consequence, the screening length is strongly influenced by the setup, i.e. whether or not semi-core states are included. In many cases, semi-core states may be considered not to contribute much to screening, and the screened exchange may overestimate the extent of screening, if the screening length is determined automatically. The alternative is to provide a

specified value: if this option is selected an additional entry field for the Thomas Fermi screening length becomes available to enter appropriate data in units of 1/\({\mathring{\mathrm{A}}}\).

For any of the functionals based on non-local Hartree-Fock exchange can be chosen from the Functional pull-down menu, a frame for Technical Settings for Non-local Exchange appears (see above Figure), which allows you to customize parameters and procedures involved in non-local exchange calculations. The procedures and options provided by this frame is discussed below:

Calculations involving non-local exchange tend to be more difficult to converge than usual DFT calculations. Therefore, any hybrid functional, screened exchange or Hartree-Fock calculation should be initialized by a DFT calculation. This is implemented and automatized when running within the MedeA environment. As a standard behavior, each task applying a non-local exchange functional as defined by Type of calculation (single point, geometry optimization, molecular dynamics) is initialized from the wave functions obtained from a single point calculation for the same system and the same computational parameters but applying a density functional. Whether the initial DFT calculation is a single point energy run or a more involved simulation can be specified for each Type of calculation by the Protocol pull-down menu.

Protocol:

For Type of calculation set to Single Point the options are:

DFT Single Point + Non-local Single Point: The non-local exchange single point energy calculation is initialized from the wave functions obtained from a DFT single point energy run.

Non-local Single Point Only: The initial DFT calculation is dismissed and the non-local exchange single point energy calculation is started from scratch.

For Type of calculation set to Structure Optimization the options are:

DFT Single Point + Non-local Structure Optimization: The non-local exchange structure optimization is initialized from the wave functions obtained from a DFT single point energy run. The main interest is the structure as obtained from a non-local exchange functional.

DFT Structure Optimization + Non-local Single Point: The structure optimization is performed on the DFT level only and accurate total energies and properties are obtained from the non-local exchange functional initialized from the DFT wave functions. This protocol is motivated by the idea that geometries may be less affected by the local or semi-local approximations of DFT, and that improvements are mainly expected for the total energy and other properties such as the band gap, magnetic moments, etc.

Non-local Structure Optimization Only: The initial DFT calculation is dismissed and the structure optimization is started on the non-local exchange level from scratch.

DFT Structure Optimization + Non-local Structure Optimization: The structure is optimized both on the DFT and non-local exchange level and results for the optimized structural data can be compared after a single calculation. Of course, the wave functions of the DFT are used to start the non-local exchange calculation.

For Type of calculation set to Molecular Dynamics the options are:

DFT Single Point + Non-local Molecular Dynamics: The wave functions of a single point energy DFT calculations serve as a starting point to initialize a non-local exchange molecular dynamics run.

DFT Molecular Dynamics + Non-local Single Point: The molecular dynamics simulation is performed on the DFT level only and for the last frame accurate total energies and properties are obtained from the non-local exchange functional as initialized from the DFT wave functions.

Non-local Molecular Dynamics Only: The initial DFT calculation is dismissed and the molecular dynamics simulation is started on the non-local exchange level from scratch.

The two-step protocols discussed above are not only applied for the simulation selected for the Type of calculation but also for some of the properties chosen from the Property frame, both available in the Calculation panel. Valence, difference, and spin charge densities, total local potentials, wave functions, and electron localization functions are provided for the second protocol step only, which for the default protocols always apply the non-local exchange level of theory. For the calculation of band structures, densities of states and optical spectra (the properties requiring additional tasks), however, the two-step protocol is applied and these properties become available both on the DFT and non-local exchange level of theory. Applying the non-local exchange functional, data are written to BandStructure.data, DensityOfStates.data and OpticalSpectra.data, whereas the corresponding data calculated from DFT are written to DFT_BandStructure.data, DFT_DensityOfStates.data and DFT_OpticalSpectra.data. Each property can be viewed and analyzed from the Analysis menu in MedeA.

Starting from corresponding DFT wave functions helps, and in many cases enables, the convergence of non-local exchange calculations. However, none of the algorithms for optimization of the electronic ground state and charge density mixing proving successful and efficient for standard DFT is applicable for optimizing the ground state for non-local exchange functionals. To this end, alternative algorithms and mixing schemes are available from

Non-local exchange algorithm:

Damped molecular dynamics: This algorithm is recommended for most cases, in particular for small band gap semiconductors and metals. The performance critically depends on the chosen Time step size. For slow convergence an increase, and for divergent-like behavior a decrease of this parameter is indicated.

Preconditioned conjugate gradient: This algorithm is recommended for insulators. The best stability is usually obtained if the number of bands equals half the number of electrons (non spin polarized case). For small gap systems it is desirable and for metals it is required to use a larger value for the number of bands (see Advanced/Restart panel). For these systems the damped MD algorithm is recommended. The stability of this algorithm depends on the Initial time step size.

Normal (blocked Davidson) + Kerker: The blocked Davidson algorithm tends to be rather slow for non-local exchange, and in many cases the Pulay mixer (default for DFT) is unable to determine the proper ground-state. Therefore, a Kerker like mixing is applied and the mixing parameter needs to be decreased to a value that allows convergence. To revert to Pulay mixing leave the Mixing parameter (Kerker) field empty.

Two parameters are provided to tune the performance of non-local exchange calculations:

Maximum angular quantum number for charge augmentation: Depending on the elements in the system different settings should be applied. The smaller the maximum angular quantum number, the faster the calculation. Testing of this parameter is recommended. The options are:

lmax = 2: This choice might be applicable for s-p systems (semiconductors), but the accuracy needs to be tested carefully.

lmax = 4: This is applicable for s-p systems and systems containing transition metals.

lmax = 6: For f electron systems the recommended choice.

lmax = 8: No truncation of the angular momentum expansion. Accurate but computationally very expensive

Reduce cutoff (FFT grid) for non-local exchange: If checked the smallest possible Fast Fourier Transformation (FFT) grid is applied for the non-local exchange part. This accelerates the calculations by roughly a factor two to three, but causes slight changes in the total energies and a small noise in the calculated forces.

2.10.7. The SCF Panel

Constructing the electronic density for a given arrangement of atoms and evaluation of the total energy, forces and stress tensors involve a self-consistent-field (SCF) calculation solving the Kohn-Sham equations. The settings needed for this step are controlled by the SCF panel and are applied for all simulations specified by the Type of calculation pull-down menu (and those properties not involving separate tasks (see section 2. Properties).

Three frames are provided: k-mesh definition, k-space integration scheme, and SCF control parameters. In MedeA VASP 5 the SCF panel is different depending on whether density functional or any of the non-local exchange functionals are selected by the Functional pull-down menu because of different requirements for the k-mesh definitions. For density functional calculations the k-mesh is specified by the options provided in the frame K-mesh in Brillouin Zone for SCF. For non-local exchange, however, two different k-meshes need to be specified: one for the evaluation of the Hartree-Fock exchange and a second one for evaluation of total energies, forces, stress tensors etc. (SCF mesh). The mesh for Hartree-Fock exchange should be chosen as small as possible, because it largely dominates the CPU time and memory demands. The SCF mesh cannot be chosen independently in VASP 5: it can only apply mesh subdivisions being an integer multiple of the subdivisions used for Hartree-Fock exchange. Therefore, the SCF mesh is derived by applying multiplication factors to the subdivisions of the Hartree-Fock mesh. The SCF panel for non-local exchange calculations exhibits two frames for k-mesh specification: the frame K-mesh in Brillouin Zone for Non-local Exchange with the usual options for the k-mesh definition and a second frame for K-mesh in Brillouin Zone for SCF for applying the multiplication factors. The two different appearances of the SCF panel are shown by the below.

The SCF panel for density functional calculations:

The SCF panel for non-local exchange calculations:

The below descriptions of parameters and options for the k-mesh definition and k-space integration are generally applicable, they reappear for instance for the density of states, optical spectra, response tensors in the DOS/Optic/Tensors panel, see section 7. The DOS/Optic/Tensors panel).

2.10.7.1. K-mesh in Brillouin Zone

There are two different input modes for specifying the k-mesh in the Brillouin zone:

Input mode:

set spacing between k-points - Uses a regular mesh of k-points in reciprocal space as equidistant as possible, defined by the Spacing of k-points entry in units of 1/\({\mathring{\mathrm{A}}}\)

set mesh parameters explicitly - Allows definition of explicit values for the k-mesh subdivisions in each of the X, Y and Z directions of k-space

Note

The k-mesh is one of the most critical parameters for an electronic structure calculation. Since convergence behavior is heavily dependent on the system, it is strongly recommended to test your models for the convergence of the properties of interest with respect to the k-mesh density.

For example, the convergence behavior of the total energy can vary between odd and even sized grids. Carefully check or use either odd or even meshes exclusively. Use a fixed k-spacing to compare calculations for models with different shapes/volumes, keeping an eye on the information provided by Actual mesh and k-spacing (see below).

Shift origin to Gamma: Shifts the origin of the k-mesh to the G-point

Use odd sized grids: Always have an uneven number of grid points in each direction

Both these options generate k-meshes with their origin at the G-point, either by shifting k-meshes with even number of mesh points in one or more directions into G or by restricting mesh parameters to odd numbers.

Actual mesh and k-spacing: This frame shows the number of k-mesh subdivisions and the actual spacing between k-points in each direction as calculated from the provided input. This frame also indicates constraints on the k-mesh due to symmetry, e.g. for cubic crystal class x=y=z, for tetragonal and hexagonal crystal classes x=y. Due to the construction of the k-meshes also less obvious constraints may be imposed (x=y=z for centered tetragonal lattices).

Note

For surface calculations, you may ensure that only 1 (or at maximum 2) k-point subdivision in the reciprocal direction of the surface normal is used. For molecular systems use a single k-point (G-point only) as generated by only 1 subdivision in each of the directions.

For non-local exchange calculations, the k-mesh for SCF must be derived from the k-mesh for the Hartree-Fock term by application of multiplication factors X axis k-mesh factor, Y axis k-mesh factor, and Z axis k-mesh factor provided in a second frame for the k-mesh definition. There is a further Actual mesh and k-spacing frame summarizing the actual settings and constraints. (It is noted that modification of the Hartree-Fock mesh data does not update the SCF k-mesh data automatically.)

2.10.7.2. Integration Scheme

The choice of the integration scheme determines how the electronic density of states is integrated. In the limit of an infinitely fine integration mesh, all integration methods should yield identical results. In reality, however, relatively coarse meshes are used to speed up computations. The most sensitive area in k-space is the area dividing occupied from unoccupied states. For semiconductors and insulators exhibiting a band gap states are either occupied or unoccupied. However, for metals a Fermi surface exists and due to the limited number of k-points fractional occupations need to be considered, occupied and unoccupied states are direct neighbors in energy and k-space. Therefore, numerical integration of the density of states of metals needs great care and advanced algorithms. Several different schemes to achieve this are implemented in VASP:

• Methfessel-Paxton, Fermi, Gaussian: Use smearing of the electronic occupation around the Fermi energy to improve the convergence of integration results with the number of k-points. The default Smearing width is 0.2 eV for Methfessel-Paxton and 0.05 eV for Fermi and Gaussian smearing.
• Tetrahedron, Tetrahedron with Blochl corrections: Use a tetrahedra decomposition of the Brillouin zone to integrate the electronic density by linear interpolation inside the tetrahedra. A correction term suggested by Bl\({\ddot{o}}\) chl overcomes to some extent the linear approximation and improves convergence with the number of k-points.

Recommendations:

Metals: Use Methfessel-Paxton for structure optimization, in particular for large cells and supercells. If default smearing produces a large entropy term (>1meV) test with varying smearing. For precise total energy calculations for small or medium size unit cells use a fine k-mesh (0.2 1/\({\mathring{\mathrm{A}}}\) or better) and Tetrahedron with Blochl corrections.

Semiconductors/Insulators: Use Tetrahedron method or Tetrahedron with Blochl corrections. If the number of k-points is too small and unit cell sizes too large to allow for the tetrahedron method to be used, use Gaussian. In this case the Smearing width, which should be smaller than for Methfessel-Paxton. If in doubt use the automated convergence tool to test.

The Smearing width: numerical parameter is used to define the width of the smearing function for Methfessel-Paxton, Fermi, and Gaussian techniques

Smearing out the electronic density around the Fermi surface region is important for metals as the determination of the precise location of the Fermi energy is numerically difficult.

The convergence of the SCF procedure may become improved by applying a smearing method for the integration of the electronic density of states, since density fluctuations between iterations during the SCF cycle are suppressed.

Order of smearing function: The order of the Methfessel-Paxton smearing

functions to be used: linear, quadratic …

2.10.7.3. SCF Control

SCF convergence: Determines the convergence of the self-consistent-field run. Convergence is reached if both the total energy and the electronic eigenvalues (band structure energies) change less than the value given in the input field (in eV) between two subsequent iterations.

Note

The default of 1E-5 eV is sufficient for most structure relaxations. Set the SCF convergence to 1E-6 or 1E-7 if the structure optimization is required to yield more precise forces/geometries (e.g. a Phonon atomic minimization of the “undisplaced” structure)

Maximum iterations: The maximum number of iterations before stopping the SCF cycle

Minimum iteration: The minimum number of iterations to do in an SCF cycle

Initial delay: number of initial steps to do without updating the wave functions (non self-consistent)

Note

In most cases, the SCF cycle converges within less than 60 steps. When starting a structure optimization from a very unrealistic input structure, the SCF may not converge during the first few geometry steps. An initial delay of 8-12 may be required for difficult surface calculations

2.10.8. The DOS/Optic/Tensors Panel

The DOS/Optic/Tensors panel controls the precision of the k-point sampling for the calculation of the electronic density of states (DOS), optical spectra and response tensors, the integration scheme for each of these tasks and parameters for defining the projection scheme, grid, energy range, number of bands, and complex shift parameters. Calculating the DOS requires a fairly dense mesh of k-points in the Brillouin zone (k-space sampling).

In MedeA VASP 5 the DOS/Optic/Tensors panel is different depending on whether density functional or any of the non-local exchange functionals are selected by the Functional pull-down menu because of different requirements for the k-mesh definitions. For density functional calculations the k-mesh is specified by the options provided in the frame K-mesh in Brillouin Zone. For non-local exchange calculations, however, the k-mesh cannot be chosen independently in VASP 5: it can only apply mesh subdivisions being an integer multiple of the subdivisions used for Hartree-Fock exchange. Therefore, the k-mesh in this panel is derived by applying multiplication factors to the subdivisions of the Hartree-Fock mesh, as discussed for the SCF.

A MedeA DOS calculation consists of two steps (four steps if a non-local exchange functional is chosen):

• A self-consistent calculation is chosen as Type of calculation to generate a converged charge density satisfying the criteria selected in the SCF panel
• If non-local exchange is selected as Functional step 1 and 2 are executed as specified by the Protocol selected from the Functional/Potential panel.
• A restart using the converged charge density from the first step and applying a different (usually finer) k-point sampling and probably a different integration technique to calculate the density of states.
• If non-local exchange is selected as Functional the DOS calculation is repeated based on this functional making use of the wave functions of step 2.

The DOS/Optic/Tensors panel for density functional calculations:

The DOS/Optic/Tensors panel for non-local exchange calculations:

The DOS/Optic/Tensors panel offers comparable options as the SCF panel for the k-mesh definition and integration scheme, but with different defaults and additional entries to achieve higher precision:

• The spacing of k-points is somewhat decreased (0.25 1/\({\mathring{\mathrm{A}}}\)) leading to a denser k-mesh
• For non-local exchange calculations the X axis k-mesh factor, Y axis k-mesh factor, and Z axis k-mesh factor are set to 2, whereas the default for SCF is 1. This leads to a denser k-mesh in this case, too.
• K-mesh for response tensors, NMR and phonons: This choice allows selecting the k-mesh for linear response tensors, NMR chemical shifts, and zone center phonons from the two available meshes. These are the meshes:

as for DOS and optics - the usually finer mesh as defined in the DOS/Optic/Tensors panel. This option is the default because linear response tensors and NMR chemical shifts are found to be quite sensitive to k-mesh size.

as for SCF - the usually coarser mesh as defined in the SCF panel

• By default, the Tetrahedron method is used as an integration scheme for the density of states and optical spectra. The Tetrahedron method with Blochl corrections is not provided, because there is no difference to the standard tetrahedron method for the DOS. The tetrahedron method is recommended for the DOS of materials with small or medium unit cell size, whereas for surfaces or large unit cells using only a few k-points the other integration techniques are more appropriate. For the calculation of optical spectra, the tetrahedron method may cause errors. Gaussian smearing or one of the integration techniques other than tetrahedron method are recommended for optical spectra, therefore. The density of states could be run separately with the tetrahedron method, if needed.

  By default, the Tetrahedron method with Blochl corrections is used as an integration scheme for response tensors and NMR chemical shifts. The Methfessel-Paxton method is not provided, because linear response is applicable for semiconductors and insulators only, and this smearing technique was found to cause substantial numerical errors.

• In contrast, for the zone center phonons the Methfessel-Paxton method is expected to yield best convergence behavior for metals and is chosen as default, therefore. For insulators and semiconductors, the Gaussian smearing method with a reduced Smearing width is recommended.

In the Parameters for DOS and Optical Spectra frame, the following options can be set to customize density of state and optical spectra simulations:

• DOS projection onto: The projection scheme applied to obtain site and angular momentum (s, p, d, f) projected density of states. For serial job execution all projection schemes (onto spherical harmonics, PAW spheres or Bader volumes, which requires an extra task to calculate the total charge density) are applicable within any functional. For parallel job execution, however, there are limitations for the choice of the projection scheme, depending on the applied functional:
• Projection onto spherical harmonics and Bader volumes requires parallelization over plane wave coefficients only (NPAR=1), the first option automatically switched on in the density of states part. Parallel calculations applying non-local exchange functionals are implemented for parallelization over bands only (NPAR=number of nodes). In this case, the parallelism of the density of states part requires the projection onto PAW spheres, which is automatically switched on. The density of states for non-local exchange functionals cannot be run in parallel projecting onto spherical harmonics. The “automatic choice” applies projection onto spherical harmonics for all density functional based DOS calculations (setting NPAR=1) and applies projection onto PAW spheres for all non-local exchange functionals (NPAR=number of nodes). The same choice is made even if calculations are executed serially.
• In summary, feasible options are:
• DFT DOS: projection onto spherical harmonics serial or parallel NPAR=1 projection onto PAW spheres serial or any parallel (default NPAR=number of nodes
• Non-local DOS: projection onto spherical harmonics serial projection onto PAW spheres serial or parallel NPAR=number of nodes
• Automatic choice:
  • DFT DOS: projection onto spherical harmonics serial or parallel NPAR=1
  • Non-local DOS: projection onto PAW spheres serial or parallel NPAR=number of nodes

Warning

projection onto PAW spheres may give zero d or f components, which actually should be non-zero. Example: a zero d component of the DOS is obtained for silicon.

• Number of grid points: This specifies the number of energy grid points for calculating and visualizing the electronic density of states and for optical spectra calculations the number of frequency points, in addition. The default is 3000 grid points, which is sufficient for most cases.
• Minimum energy and Maximum energy (in eV) are defined relative to the Fermi level and define an energy “window”. This option may be used to cut off low-lying peaks of semi-core states, which may decrease the resolution of the energy mesh for valence electron states.
• Number of bands: The number of bands used for the calculation of the density of states and optical spectra is an important parameter. In particular for optical spectra a rather large number of unoccupied bands is required. If the entry field is empty the default number of bands as shown above the entry field is applied, which is a much larger number for optical spectra than for the density of states only. If calculations are run parallel, this number may become slightly increased such that the number of bands becomes a multiple of the number of processors working in parallel.
• Shift parameter: The small complex shift used in the Kramers-Kronig transformation smoothing the real part of the dielectric function. The default value is appropriate for most cases. Only for systems with very small band gap (about two times the specified shift parameter) the static value of the dielectric function may become inaccurate. In such cases the shift parameter should be decreased and the number of grid points should be raised to values of about 2000.
• Print optical matrix elements: Provides a listing of optical matrix elements (oscillator strengths) for all transitions in a separate output file, allowing detailed analysis.

The Parameters for NMR chemical shifts frame provides access to

• Planewave cutoff for NMR: to specify an increased planewave cutoff in eV required for NMR chemical shift calculations. The default value shown above is increased by 50% above the usual default cutoff.

The Parameters for Zone Center Phonons frame holds further input parameters for this property calculation

• Displacement: The magnitude of displacements of atoms from their equilibrium positions in Angstrom units for calculating finite differences of forces.
• Number of displacements: The number of central difference displacement values. Multiples of the above displacement value are applied in a positive and negative direction. Using more than 1 displacement may increase the accuracy. The value of 0 would displace only in one direction, which is not recommended!

The Parameters for Total Charge Density frame allows you to

• Refine the Fourier grid by: a percentage by which to increase the Fourier grid above the default value for the total charge density and Bader decomposition calculations. The default Fourier grid is the minimum one avoiding aliasing errors (corresponding to accurate precision).

2.10.9. The Band Structure Panel

A band structure is a plot of the electronic eigenvalues as a function of the electronic momentum k. It is also referred to as the electronic dispersion relation. In principle, the band structure is a 3 dimensional scalar field e_i(k), where the vector k takes on all values within the first Brillouin zone of the crystal and i labels the bands. In practice, the symmetry of a crystal reduces the number of non-equivalent k-vectors (or k-points) considerably. One usually plots the band energies along a path connecting points of high symmetry within the irreducible sector of the Brillouin zone of the crystal.

MedeA provides standard paths for all lattice types. Making use of standard paths only two input parameters need to be set, which are

Maximum number of points: The total number of k-points calculated by VASP for plotting the band structure, aligned along the default path as equidistant as possible.

Number of k-points per task: Use this option to limit the number of k-points per VASP task. Main memory requirement and compute time depend on the number of k-points used. MedeA splits a band structure job into separate tasks with the Number of k-points per task k-points in each task.

The standard paths cover most of the high symmetry points and directions, but not necessarily all of them. In addition, MedeA :sup:`` offers a way to set up a user-defined path for displaying the band structure by modification of the standard path. To access this feature you may click

Manually define path: This displays the standard path of k-points available for the given system, i.e. the high-symmetry points (vertices) of the path and the number of points on the line segments between them. The path can be modified by editing the list of k-point vertices and by changing the number of line segment points to be calculated and displayed.

Vertex: A high symmetry point inside the Brillouin zone of the given crystal. The band structure displays the electronic energy as a function of the electronic momentum k on a sequence of lines connecting points (vertices) of high symmetry in the Brillouin zone. Each line segment represents a different direction in k-space.

Note

Axis definitions can vary between publications, resulting in varying coordinates for symmetry points. If in doubt, check out MedeA’s chapter on the Brillouin zone definition.

Line segments: The number of points to calculate and display on a line between two vertexes

To modify and edit Vertex and Line segments a number of operations can be used:

• Click on a vertex point and select a different one from a list of vertex points alcreated defined to replace it. The list of vertex points becomes available in the popup-menu upon clicking.
• Click on a vertex point and select –Delete this point– in order to delete it from the path
• Click on a vertex point and select –Define a new point– in order to replace the vertex point by a new one not yet in the list of vertex points. A panel opens up to enter a Label and three coordinates X, Y and Z.
• In order to add a new vertex point at the bottom of the path, click on –Add a new point– and either chose one of the points in the list or define a new one by –Define a new point–.
• Modify the number of points in each line segment accordingly
• The standard path can be recovered at any time by clicking Reset to default path.
• If a non-local exchange functional is selected by the Functional menu, an additional choice appears at the top of the Band structure panel:
• K-mesh for underlying SCF: This defines which k-mesh to apply for calculating the self-consistent charge density. For non-local exchange functionals the k-mesh applied for non-local exchange and for SCF must be identical for calculating the band structure. If the X, Y, and Z axis k-mesh factors specified in the SCF panel is not set to 1, this choice allows you to decide whether the k-mesh for non-local exchange or that one for SCF are applied as a basis for the band structure evaluation. It is noted that large k-meshes for the non-local exchange increase drastically the computational demands.

2.10.10. The Advanced/Restart Panel

The Advanced/Restart panel summarizes technical control settings, restart settings and options for including strong correlations in the crystal Hamiltonian. For standard calculations, it is recommended to keep these defaults. A detailed description of the available settings is given below:

2.10.10.1. Technical Settings

Algorithm:

Normal (blocked Davidson) - Recommended default

Fast - Starts with blocked Davidson and switches to RMM-DIIS after a number of steps

Very fast: The RMM-DIIS algorithm reduces the number of normalization steps considerably and is therefore much faster than the Davidson algorithm for large systems and on workstations with a small memory bandwidth.

Damped molecular dynamics: This algorithm is a powerful alternative for difficult cases where Normal (blocked Davidson) tends to fail. The performance critically depends on the chosen Time step size. For slow convergence an increase, and for divergent-like behavior, a decrease of this parameter is indicated.

Preconditioned conjugate gradient: The best stability is usually obtained if the number of bands equals half the number of electrons (non-spin polarized case). The stability of this algorithm depends on the Initial time step size.

Normal (blocked Davidson) + Kerker: For cases where the Pulay mixer (applied for Normal, Fast, and Very fast) is unable to determine the proper ground state. A Kerker like mixing is applied and the mixing parameter needs to be decreased to a value that allows convergence. To recourse to Pulay mixing leave empty the entry field for Mixing parameter (Kerker).

Note

RMM-DIIS stands for Residual Minimization Method with Direct Inversion of the Iterative Subspace.

It does not have as large a radius of convergence as the Davidson method, so it may have convergence problems for some systems, particularly if starting from guessed wave functions.

We recommend testing the ‘Fast’ or ‘Very Fast’ RMM-DIIS algorithms together with real space projection for large systems.

Extra fine augmentation grid: Check this option to use an extra fine augmentation grid to yield very accurate forces. This option is recommended for eliminating numerical noise in complex Phonon calculations. Higher accuracy is achieved by defining an additional super fine Fourier grid doubling the sampling of the fine grid used to evaluate the augmentation charges.

Do not use symmetry: VASP does not make use of symmetry, if this option is checked. Not using symmetry increases dramatically the computational demands, in particular because the k-point mesh cannot be restricted to a symmetry-irreducible set. This option should be considered only for very specific cases, such as spin-orbit relativistic (e.g. for magnetic anisotropies), non-collinear magnetic calculations or DFT + U calculations. Note: Molecular dynamics calculations are always run without symmetry, independent of this setting.

Use ScaLAPACK: VASP can make use of the ScaLAPACK library (Scalable Linear Algebra PACKage) providing high-performance linear algebra routines for parallel distributed memory machines. Currently, MedeA VASP does not use ScaLAPACK per default. Use of ScaLAPACK is particularly recommended for GW and ACFDT-RPA calculations.

Number of bands: The number of bands to include in the calculations. At least one empty band should be present. VASP issues a warning in if this is not the case. The default number of bands is displayed above the entry field.

To check this parameter, perform several calculations for a fixed potential (ICHARG=12 in Add to Input) with an increasing number of bands. An agreement of the total energies respective of 1e ^-6 should be obtained in 10-15 iterations

Mind that the RMM-DIIS scheme is more sensitive to the number of bands then the default CG algorithm.

The actual value is defined by the keyword NBANDS and can be found in OUTCAR.out.

2.10.10.2. Magnetism

Spin interpolation:

Vosko-Wilk-Nusair - consistent with both LDA and GGA calculations

Barth-Hedin - consistent with LDA calculations

The earlier Bart-Hedin algorithm was used to set up the exchange correlation potential in LDA. With the arrival of the GGA formalism, the Vosko-Wilk-Nusair implementation has become a standard.

Total magnetic moment: Allows to constrain the total magnetic moment to the specified value in mB during a self-consistent spin-polarized calculation.

Spin quantization axis: This becomes available only if Spin-orbit magnetic is chosen for Magnetism. Magnetic anisotropies can be calculated by varying the spin quantization axis while keeping the non-collinear magnetic moments of the atoms constant.

2.10.10.3. Element/site specific choices

The pulldown menu

Enable choices specific for - enables to switch between elements and sites, which triggers the behavior in several parts of the GUI. In the Functional/Potential panel either one potential can be chosen for each element, or different potentials can be chosen for different sites occupied by atoms of the same element. This could be used, for instance, to saturate dangling bonds with different types of hydrogen atoms.

Another application is in treating strong correlation by the LDA + U approach (see section iii below): In elements specific mode one set of U and J parameter can be chosen for one element, even if the atoms are in different chemical environments. In sites specific mode U and J parameters can be set for each site, in the example below different parameters for Fe²⁺ and Fe³⁺ in Fe ₃O ₄.

File return: Enables customization of the extent to which files are returned from the TaskServer back to the job directory on the JobServer. The options are:

Normal - returns the main output files only, avoiding very large files

Normal + vasprun.xml - same as Normal but returns in addition the vasprun.xml file, which can be quite large

All files - returns all files created by a VASP run

2.10.10.4. Initial Conditions and Restart

Initial wave functions:

From scratch - reinitializes the wave functions at the beginning of a SCF run

Read in from previous run - use wave functions from a previous calculation.

Choosing From scratch, there are two options to proceed with initialization:

Initialize wave functions with:

Random number - uses random numbers for initialization

Jellium - the Jellium wave function fills the plane waves with the lowest kinetic energy for a constant potential

Choosing Read in from previous run there are two possible ways how the new wave functions setup is initialized from the restart wave functions:

Restart calculation with:

constant energy cutoff - If the size and shape of the cell or the energy cutoff has changed with respect to the previous run, this option redefines the set of plane waves according to the new cell. This option is recommended for any total energy based cell size and shape optimizations and also for convergence tests.

constant basis set - This should only be applied if an optimization of cell volume and shape is to be continued with a consistent basis set. The continuation will then be equivalent to optimizing in a single calculation. The setup will not be adapted to the shape of the cell.

Note

Note that for a fixed energy cutoff the number and setup of plane waves at a given k-point depends implicitly on the cell parameters, which is addressed by the above settings.

Initial charge density:

From initial wave functions: If the initial wave functions are initialized from scratch, it is recommended to use a superposition of atomic charge densities

Read in from previous run: If suitable wave functions and/or charge densities are available, using them to initialize the SCF run will speed up convergence considerably

Atomic charge densities: The initial charge density is constructed as a superposition of atomic charge densities

Fix the charge density fixes charge density to its initial value for the whole of the SCF calculation. Fixing the charge density can be useful for calculating the band structure, DOS, or for running a Harris-Foulkes functional calculation (using the atomic densities)

2.10.10.5. Strong Correlation

Standard LDA or GGA - Standard density functional approach without L(S)DA+U

Simplified LSDA+U - The approach suggested by Dudarev, Savrasov, Humphreys & Sutton [62]: only one parameter U-J is required to define the on-site Coulomb interaction

Rotationally invariant LSDA+U The approach suggested by Liechtenstein & Zaanen [63] : both parameters U and J need to be specified

Rotationally invariant LDA+U - Like the LSDA+U variant above, however, with a different definition of the so-called double counting energy

Note

The unscreened electron-electron interaction can be written in terms of Slater integrals. Application of Slater integrals as calculated from atomic wave functions results in a large overestimation of the true interaction, since the Coulomb interaction in condensed systems is screened. Therefore, in practice these integrals are evaluated in terms of the Coulomb and exchange parameters U and J.

These parameters are adjusted to obtain agreement with experimental data such as equilibrium volume, magnetic moments, band gap, or crystal structure. This approach is known as L(S)DA+U method. Despite its name, it can equally well be applied within GGA.

2.10.11. The Add to Input Panel

The MedeA VASP graphical user interface does not expose all possible VASP parameters. To give access to all the features of VASP, which are not explicitly part of the interface, VASP keywords and parameters can be entered as additional lines of the INCAR file. This is accomplished by providing the Add to Input panel

Additional Input Lines

NBANDS = 250

As the VASP interface remembers settings for the next calculation, the added input lines might not make sense during the next run. To give you a reminder, the panel is highlighted as Add to Input in red, whenever Additional Input Lines have been specified:

Note

VASP Keywords and parameters are described in the The VASP Manual` Note that VASP uses the first instance of a keyword, thus settings defined by the user in Add to Input will override settings defined elsewhere in the VASP interface panels.

The modified input files can be inspected in the Preview panel.

2.10.12. The Preview Input Panel

File INCAR

# SCF input for VASP
# Note that VASP uses the FIRST occurrence of a keyword
SYSTEM = (Al)4 (Fm-3m) ~ Al (VASP)
PREC = Accurate
ENCUT = 500.000
IBRION = 2
NSW = 100
ISIF = 3
ALGO = Normal (blocked Davidson)
...

Copy to clipboard

This panel lets you preview VASP input files for the simulation specified by Type of calculation before actually launching a job. Select the type of file you would like to preview from the Files menu:

• POTCAR - lists the potentials used for each atom in the system
• script - job submission script
• INCAR - the main VASP input file
• KPOINTS - contains k-mesh information
• POSCAR - contains atomic positions

This feature is often used to preview and check settings added by hand using the Add to Input capability. Note that only the first occurrence of a given keyword will be used and options set elsewhere in the graphical interface are overwritten by Add to Input.

[1]	Georg Kresse and Jürgen Furthmüller, “Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set,” Computational Materials Science 6 (1996): 15.

[2]	Walter Kohn and L J Sham, “Self Consistent Equations Including Exchange and Correlation Effects,” Physical Review A 140, no. 4 (1965): 1133-1138.

[3]	Georg Kresse and D Joubert, “From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method,” Physical Review B 59, no. 3 (1999): 1758.

[4]

Georg Kresse and Jürgen Furthmüller, “Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set,” Physical Review B 54, no. 16 (1996): 11169.G Kresse and J Furthmüller, “Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set,” Computational Materials Science 6 (1996): 15.

[5]	Georg Kresse and D Joubert, “From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method,” Physical Review B 59, no. 3 (1999): 1758.

[6]	Peter E Bl\({\ddot{o}}\) chl, “Projector Augmented-Wave Method,” Physical Review B 50, no. 24 (December 1994): 17953-17979.

[7]	S Nosé, “Constant Temperature Molecular Dynamics Methods,” Progress of Theoretical Physics Supplement 103 (1991): 1.

[8]	(1, 2) M. Parrinello, A. Rahman, “Crystal Structure and Pair Potentials: A Molecular-Dynamics Study”, Physical Review Letter 45, (1980): 1196

[9]	(1, 2) M. Parrinello, A. Rahman, “Polymorphic transitions in single crystals: A new molecular dynamics method”, Journal of Applied Physics 52, (1981): 7182

[10]	J. Paier, M. Marsman, and G. Kresse, “Dielectric properties and excitons for extended systems from hybrid functionals”, Physical Review B 78, (2008): 121201 (R)

[11]	L. Hedin, “New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem”, Physical Review 139, (1965): A796

[12]	M.S. Hybertsen, S.G. Louie, “Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies”, Physical Review B 34, (1986): 5390

[13]	T. Sander, E. Maggio, G. Kresse, “Beyond the Tamm-Dancoff approximation for extended systems using exact diagonalization”, Physical Review B 92, (2015): 045209

[14]	M. Shishkin, G. Kresse, “Implementation and performance of the frequency-dependent GW method within the PAW framework”, Physical Review B 74, (2006): 035101

[15]	M. Shishkin, G. Kresse, “Self-consistent GW calculations for semiconductors and insulators”, Physical Review B 75, (2007): 235102

[16]	M. Shishkin, M. Marsman, G. Kresse, “Accurate Quasiparticle Spectra from Self-Consistent GW Calculations with Vertex Corrections”, Physical Review Letters 99, (2007): 246403

[17]	J. Harl, G. Kresse, “Accurate Bulk Properties from Approximate Many-Body Techniques”, Physical Review Letter 103, (2009): 056401

[18]	J. Harl, L. Schimka, G. Kresse, “Assessing the quality of the random phase approximation for lattice constants and atomization energies of solids”, Physical Review B 81, (2010): 115126

[19]	Pekka Pyykk\({\ddot{o}}\) , “Year-2008 nuclear quadrupole moments”, Mol. Phys. 106, (2008): 1965

[20]

W. Tang, E. Sanville, and G. Henkelman, “A grid-based Bader analysis algorithm without lattice bias”, Journal of Physics: Condensed Matter 21, (2009): 084204. E. Sanville, S. D. Kenny, R. Smith, and G. Henkelman, “An improved grid-based algorithm for Bader charge allocation”, Journal of Computational Chemistry 28, (2007): 899-908. G. Henkelman, A. Arnaldsson, and H. J\({\acute{o}}\) nsson, “A fast and robust algorithm for Bader decomposition of charge density”, Computational Materials Science 36, (2006): 254-360

[21]	C.J. Pickard, F. Mauri, “All-electron magnetic response with pseudopotentials: NMR chemical shifts”, Physical Review B 63, (2001): 245101

[22]	J.R. Yates, C.J. Pickard, F. Mauri, “Calculation of NMR chemical shifts for extended systems using ultrasoft pseudopotentials”, Physical Review B 76, (2007): 024401

[23]	K. Mathew, R. Sundararaman, K. Letchworth-Weaver, T.A. Arias, R.G. Hennig, “Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways”, Journal of Chemical Physics 140, (2014): 084106

[24]	John P Perdew and Alex Zunger, “Self-Interaction Correction to Density-Functional Approximations for Many-Electron Systems,” Physical Review B 23, no. 10 (1981): 5048-5079.

[25]	D M Ceperley, “Ground State of the Electron Gas by a Stochastic Method,” Physical Review Letters 45, no. 7 (August 1980): 566-569.

[26]	R Armiento and AE Mattsson, “Functional Designed to Include Surface Effects in Self-Consistent Density Functional Theory,” Physical Review B 72, no. 8 (2005): 085108.

[27]

John P Perdew, A Ruzsinszky, GI Csonka, OA Vydrov, GE Scuseria, et al., “Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces,” Physical Review Letters 100, no. 13 (2008): 136406. John P Perdew et al., “Perdew Et Al. Reply,” Physical Review Letters 101 (2008): 239702. AE Mattsson, R Armiento, and TR Mattsson, “Comment on ‘Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces’,” Physical Review Letters 101, no. 23 (2008): 239701. John P Perdew et al., “Erratum: Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces [Phys. Rev. Lett. 100, 136406 (2008)],” Physical Review Letters 102, no. 3 (2009): 39902.

[28]	John P Perdew, Kieron Burke, and Matthias Ernzerhof, “Generalized Gradient Approximation Made Simple,” Physical Review Letters 77, no. 18 (October 1996): 3865-3868.

[29]	B Hammer, L B Hansen, and Jens K Norskov, “Improved Adsorption Energetics Within Density-Functional Theory Using Revised Perdew-Burke-Ernzerhof Functionals,” Physical Review B 59, no. 11 (1999): 7413-7421.

[30]	A D Becke, “Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior,” Physical Review A 38, no. 6 (1988): 3098.

[31]

J. Klimes, D. R. Bowler, and A. Michaelides, “Van der Waals density functionals applied to solids”, Physical Review B 83, (2011): 195131; J. Klimes, D. R. Bowler, and A. Michaelides, “Chemical accuracy for the van der Waals density functional”, Journal of Physics: Condensed Matter 22, (2010): 022201.

[32]	J. Wellendorff, K.T. Lundgaard, A. Mogelhoj, V. Petzold, D.D. Landis, J.K. Norskov, T. Bligaard, K.W. Jacobsen, “Density functionals for surface science: Exchange-correlation model development with Bayesian error estimation error estimation”, Physical Review B 85, (2012): 235149

[33]	I. Hamada, “van der Waals density functional made accurate”, Physical Review B 89, (2014): 121103

[34]	K. Lee, E. D. Murray, L. Kong, B. I. Lundqvist, and D. C. Langreth, “Higher-accuracy van der Waals density functional”, Physical Review B 82, (2010): 081101

[35]	M. Dion, H. Rydberg, E. Schr\({\ddot{o}}\) der, D. C. Langreth, and B. I. Lundqvist, “Van der Waals Density Functional for General Geometries”, Physical Review Letters 92, no. 24 (2004): 246401-1

[36]	H. Peng, Z-H. Yang, J.P. Perdew, J. Sun, “Versatile van der Waals Density Functional Based on a Meta-Generalized Gradient Approximation”, Physical Review X 6, (2016): 041005

[37]	O.A. Vydrov, T. Van Voorhis, “Nonlocal van der Waals Density Functional Made Simple”, Physical Review Letter 103, (2009): 063004

[38]	R. Sabatini, T. Gorni, S. de Gironcoli, “Nonlocal van der Waals density functional made simple and efficient”, Physical Review B 87, (2013): 041108(R)

[39]	J. P. Perdew, A. Ruzsinszky, G. I. Csonka, L. A. Constantin and J. Sun, “Workhorse Semilocal Density Functional for Condensed Matter Physics and Quantum Chemistry”, Physical Review Letters 103, (2009): 026403.

[40]	J. Sun, M. Marsman, G. Csonka, A. Ruzsinszky, P. Hao, Y.-S. Kim., G. Kresse, and J. P. Perdew, “Self-consistent meta-generalized gradient approximation within the projector-augmented-wave method”, Physical Review B 84, (2011): 035117.

[41]	J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, “Climbing the Density Functional Ladder: Nonempirical Meta-Generalized Gradient Approximation Designed for Molecules and Solids”, Physical Review Letters 91, (2003): 146401

[42]	J. Sun, R.C. Remsing, Y. Zhang, Z. Sun, A. Ruzsinszky, H. Peng, Z. Yang, A. , U. Waghmare, X. Wu, M.L. Klein, J.P. Perdew, “Accurate first-principles structures and energies of diversely bonded systems from an efficient density functional”, Nature Chemistry 8, (2016): 831

[43]	J. Sun, B. Xiao, A. Ruzsinszky, “Communication: Effect of the orbital-overlap dependence in the meta generalized gradient approximation”, Journal of Chemical Physics 137, (2012): 051101

[44]	J. Sun, R. Haunschild, B. Xiao, I.W. Bulik, G.E. Scuseria, J.P. Perdew, “Semilocal and hybrid meta-generalized gradient approximations based on the understanding of the kinetic-energy-density dependence”, Journal of Chemical Physics 138, (2013): 044113

[45]	Y. Zhao and D. G. Truhlar, “A new local density functional for main-group thermochemistry, transition metal bonding, thermochemical kinetics, and noncovalent interactions”, Journal of Chemical Physics 125, (2006): 194101.

[46]	A. D. Becke and E. R. Johnson, “A simple effective potential for exchange”, Journal of Chemical Physics 124, (2006): 221101; F. Tran and P. Blaha, “Accurate Band Gaps of Semiconductors and Insulators with a Semilocal Exchange-Correlation Potential”, Physical Review Letters 102, (2009): 226401.

[47]

(1, 2) Jochen Heyd, Gustavo E Scuseria, and Matthias Ernzerhof, “Hybrid Functionals Based on a Screened Coulomb Potential,” Journal of Chemical Physics 118, no. 18 (2003): 8207. J Heyd, GE Scuseria, and M Ernzerhof, “Erratum: “Hybrid Functionals Based on a Screened Coulomb Potential, [J. Chem. Phys. 118, 8207 (2003)],” Journal of Chemical Physics 124, no. 21 (2006): 9906.

[48]	Matthias Ernzerhof and Gustavo Scuseria, “Assessment of the Perdew-Burke-Ernzerhof Exchange-Correlation Functional,” Journal of Chemical Physics 110, no. 11 (1999): 5029-5036.

[49]	Carlo Adamo and Vincenzo Barone, “Toward Reliable Density Functional Methods Without Adjustable Parameters: the PBE0 Model,” Journal of Chemical Physics 110, no. 13 (April 3, 1999): 6158.

[50]	A D Becke, “Density-Functional Exchange-Energy Approximation with Correct Asymptotic Behavior,” Physical Review A 38, no. 6 (1988): 3098.

[51]	Chengteh Lee, Weitao Yang, and Robert G Parr, “Development of the Colle-Salvetti Correlation-Energy Formula Into a Functional of the Electron Density,” Physical Review B 37, no. 2 (January 1988): 785-789.

[52]	S H Vosko, L Wilk, and M Nusair, “Accurate Spin-Dependent Electron Liquid Correlation Energies for Local Spin Density Calculations: a Critical Analysis,” Canadian Journal of Physics 58, no. 8 (August 1980): 1200-1211.

[53]	P J Stephens, F J Devlin, C F Chabalowski, and M J Frisch, “Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields,” Journal of Physical Chemistry 98, no. 45 (November 1994): 11623-11627.

[54]	S. Grimme, J. Antony, S. Ehrlich, and S. Krieg, “A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements H-Pu”, Journal of Chemical Physics 132, (2010): 154104-1-19.

[55]	S. Grimme, S. Ehrlich, and L. Goerigk, “Effect of the damping function in dispersion corrected density functional theory”, Journal of Computational Chemistry 32, (2011): 1456-1465.

[56]	A. Tkatchenko and M. Scheffler, “Accurate Molecular Van DerWaals Interactions from Ground-State Electron Density and Free-Atom Reference Data”, Physical Review Letters 102, (2009): 073005.

[57]	A. Tkatchenko, R. A. Di Stasio, R. Car, and M. Scheffler, “Accurate and Efficient Method for Many-Body van der Waals Interactions”, Physical Review Letters 108, (2012): 236402.

[58]	A. Ambrosetti, A.M. Reilly, R.A. DiStasio, “Long-range correlation energy calculated from coupled atomic response functions”, Journal of Chemical Physics 140, (2014): 018A508

[59]	T. Bucko, S. Leb\({\grave{e}}\) gue, T. Gould, J.G. \({\acute{a}}\) ngy\({\acute{a}}\) n, “Many-body dispersion corrections for periodic systems: an efficient reciprocal space implementation”, Journal of Physics: Condensed Matter 28, (2016): 045201

[60]	S.N. Steinman, C. Corminboeuf, “A generalized-gradient approximation exchange hole model for dispersion coefficients”, Journal of Chemical Physics 134, (2011): 044117

[61]	S. Grimme, “Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction”, Journal of Computational Chemistry 27, (2006): 1787-1799.

[62]	S L Dudarev, S Y Savrasov, C J Humphreys, and A P Sutton, “Electron-Energy-Loss Spectra and the Structural Stability of Nickel Oxide: an LSDA+U Study,” Physical Review B 57, no. 3 (January 1998): 1505-1509.

[63]	A I Liechtenstein and J Zaanen, “Density-Functional Theory and Strong Interactions: Orbital Ordering in Mott-Hubbard Insulators,” Physical Review B 52, no. 8 (August 1995): R5467-R5470.

download:pdf