analyze_dft

This programs is able to analyze/evaluate and prepare/manage several different kinds of VASP DFT single point property calculations. Currently supported are: Bader partial charges, partial densities of states, STM image generation and core level energies/shifts. Since each of these options offers several possible keywords to modify the analysis, this program as a two-layer structure. First, only the calculation type is chosen. Secondly, the exact options for that calculation type are chosen.

The options are (with links to the sections below):

• -bader Evaluates a Bader charge calculation.
• -ddec6 Evaluates a DDEC6 charge calculation.
• -pdos Evaluates a partial density of states calculation.
• -stm Performs the generation of STM images for surface slabs.
• -cls Manages the setup and evaluation of core level (shift) calculations for all atoms of an element.

Bader partial charges

This option evaulates Bader charge calculations and directly gives the partial charge per atom in e and files charges.pdb (to visualize the charges in vmd or similar programs) and POSCAR_charge better comparison of structure and partial charges.

Before applying this program, perform a Bader charge calculation as described in the Bader charge wiki page.

If everything is completed (including application of the chgsum.pl script and the bader program!), the file ACF.dat should have been generated in the calculation folder. This file, however, only shows the total number of valence electrons assigned to each atom and not the actual partial charge!

Now, simply apply the program by typing (ACF.dat and POSCAR need to be present in the folder):

analyze_dft -bader

and the analysis is done in usually less than a second. In the command line output, there is first given an information about the number of valence electrons per element, e.g.:

 Number of elements in POSCAR:      2
  Used number of valence electrons for the included elements:
   -Ni:   10.000000
   -H :    1.000000

Here, for each element included into the POSCAR, the number of valence electrons of the default POTCAR file as recommended by the VASP developers on their VASP-Wiki page is assumed and based on this the partial charge is calculated by the difference to the actual electron number in the ACF.dat file (if more electrons, the charge is negative, if less electrons, the charge is positive). If you chose another POTCAR file, for example Ni_pv with 16 valence electrons, simply add the valence keyword to the program, e.g.:

analyze_dft -bader -valence:Ni=16

This can be done for an arbitrary number of elements at once.

After evaluation, the average charges of all elements in the system are also printed out, giving hints to some general trends, e.g.:

 The resulting average charges are:
   -Ni:      0.002242
   -H :     -0.224199

Further, a number of output files is generated:

• bader_charges.dat contains the charges of all atoms in a simple data file, with one line per charge; the atom index in the first column and the respective charge in the second column.
• charges.pdb contains a pdb file of the current structure as well as the partial charges written in the column after the three coordinate columns. This file can be used to visualize the charges, e.g., with vmd. Open the pdb file and go to Graphics -> Representations -> Coloring Method -> Occupancy. Then, the atoms should be colored such that negatively charged atoms are red and positively are blue (almost neutral ones are white). The range of the colors can be adapted as well: Go to Trajectory in the Representations window and change there the values in the Color Scale Data Range.
• POSCAR_charge contains the POSCAR file of the system, with the charges as fourth column after the coordinates.

DDEC6 partial charges

Analogously to Bader charges, the result of a DDEC6 partial charge calculation with the Chargemol program (see wiki entry) can be evaluated.

For this, only the Chargemol output file DDEC6_tempered_net_atomic_charges.xyz needs to be present. Then, simply invoke the analyze_dft program:

analyze_dft -ddec6

The average charges of all element are written to the command line; further the different output files for visualization are written:

• ddec6_charges.dat contains the charges of all atoms in a simple data file, with one line per charge; the atom index in the first column and the respective charge in the second column.
• charges.pdb contains a pdb file of the current structure as well as the partial charges written in the column after the three coordinate columns. This file can be used to visualize the charges, e.g., with vmd. Open the pdb file and go to Graphics -> Representations -> Coloring Method -> Occupancy. Then, the atoms should be colored such that negatively charged atoms are red and positively are blue (almost neutral ones are white). The range of the colors can be adapted as well: Go to Trajectory in the Representations window and change there the values in the Color Scale Data Range.
• POSCAR_charge contains the POSCAR file of the system, with the charges as fourth column after the coordinates.

Partial densities of states

Fast and easy evaluation of VASP density of states (DOS) calculations can be done with this option. How a DOS calculation can be done, can be see in the respective (wiki page)[https://github.com/Trebonius91/utils4VASP/wiki/Density-of-States]. Note that the atom/orbital resolution must be turned on before with the LORBIT=11 keyword!

The POSCAR and DOSCAR files of the DOS calculation are needed for the evaluation.

Simply execute the program with one of the following keywords, determining the parts of the system that shall be evaluated:

• -element=[list of elements] Choose all atoms of one or several elements in the system which DOS shall be summed up. Example: -element=C,H or -element=all (choosing all atoms in the system).
• -index=[list of indices] Choose certain atoms based on their index (starting from 1). Example: -index=15,283,302.
• -orbital=[identifiers] After the atoms were chosen with -element or -index, choose the orbitals you want to consider in the DOS. In the simplest case, all orbitals can be chosen withe the option -orbital=all, or one or several of them can be chosen, for example -orbital=s,px,d (choose the s-orbital, the px orbital and all d orbitals of the chosen atoms). Thus, you can either choose all orbitals of a certain angular momentum (s, p or d) or make the choice more narrow (with respect to p and d orbitals): px, py, pz for p-orbitals, dxy, dyz, dxz, dz2, dx2y2 for d-orbitals.

Note that always atoms and orbitals must be chosen explicitly! Else, the program exits with an error. If you simply want to add up the DOS of the whole system, give the keywords -element=all and -orbital=all.

After execution, the partial DOS of the chosen atoms and orbitals is written to file dos_partial.dat. If the calculation was spin-polarized, the file will contain two columns: alpha and beta. Further, a second file dos_part_sum.dat with the sum of both spins is written in that case. Further, the total DOS is written to files dos_total.dat and dos_tot_sum.dat with the same logic as for the partial DOS.

STM image generation

This options generates STM pictures from a given PARCHG file of a surface slab (see in the STM section how to obtain this file).

Only the PARCHG file is needed as input for an arbitrary system. Two modes can be used: constant height, in which the unit cell is essentially cut through parallel to the x-y axis at the given z-value and constant current, in which the first z-value (looked from above) at which the partial charge density is equal or larger the given value is determined for each pixel in the image. A number of command line arguments are available for the program:

• -mode=[word] Select if the constant height or constant current mode shall be used. -mode=height activates the constant height mode, -mode=current activates the constant current mode.
• -pos=[value] For a constant height mode calculation, give the height (z-value) at which the local partial charge density shall be probed (in Angstroms).
• -dens=[value] For a constant current mode calculation, give the partial charge density to which the STM tip shall be lowered for each pixel of the STM picture. Usually, values between 0.01 and 0.1 give reasonable results (comparable to experimental STM images). The larger the value, the smaller and sharper the contour of the adsorbate or the substrate surface will be.
• -repeat_x=[number] How often the STM image of the given unit cell shall be repeated in x-direction (or more general, including nonorthogonal unit cells: the a-axis) to show the periodic pattern more clearly. Default: 1
• -repeat_y=[number] Similar to repeat_x, along the y-axis (or b-axis) of the unit cell.
• -grid_x=[number] The number of STM grid points (or pixels) that shall be calculated along the x-axis (or more general: a-axis) of the unit cell in PARCHG. If repeat_x is larger than 1, the total number of grid points of x will be multiplied by it (but not recalculated, thus not increasing the calculation effort significantly). Default: 1
• -grid_y=[number] Similar to grid_x, along the y-axis (or b-axis) of the unit cell.
• -gauss_width=[value] For smearing of the STM picture, to imitate an experimental STM image with limited resolution available. The given value is the width of the smearing Gaussian at half maximum (FWHM), in Angstroms. Larger values will smear the STM image out, but also increase the computation effort since more PARCHG data points need to be included in each STM grid point. Default: 0.3

The program first reads in the PARCHG file, with a value of the partial charge density at each grid point.

In the constant height mode, the x-y plane (or more general: a-b plane) at a given z-value is divided into the number of grid-points defined by grid_x and grid_y. Then, all grid points stored in PARCHG are included and weighted in a Shepard interpolation kind of interpolation. The range of indices $R_i, i=x,y,z$ included into this sum are determined from the smearing Gaussian, being of the form $\exp(-a x^2)$, by a cutoff calculate divided by the number axis length $L_i$ and the number of PARCHG points along this axis $P_i$:

$$N_i = \frac{\sqrt{\frac{\ln(100)}{a}}}{L_i \cdot P_i}$$.

Then, the current local charge density (equal to the STM intensity in line of the Tersoff-Haman methodology) is calculated by first determining the nearst PARCHG grid point $\mathbf{P_N}$ and numerically integrating over the volume determined by $N_i$ and the Gaussian weighting function and then divided by the sum of weights:

$$ S = \sum\limits_{i=P_{N,x}-N_x}^{P_{N,x}+N_x} \sum\limits_{j=P_{N,y}-N_y}^{P_{N,y}+N_y} \sum\limits_{k=P_{N,z}-N_z/2}^{P_{N,z}+N_z/2} \exp\left( -a (\mathbf{q}-\mathbf{q}_{ijk}) \right)^2 C(ijk) $$

$$ W = \sum\limits_{i=P_{N,x}-N_x}^{P_{N,x}+N_x} \sum\limits_{j=P_{N,y}-N_y}^{P_{N,y}+N_y} \sum\limits_{k=P_{N,z}-N_z/2}^{P_{N,z}+N_z/2} \exp\left( -a (\mathbf{q}-\mathbf{q}_{ijk}) \right)^2 $$

$$I_{STM,xyz} = \frac{S}{W}$$

It can be seen that the included range along z is half that of x and y, since one can assume that the resolution along z is higher than along x and y (the tip ends at a certain height and will mostly interact with the local partial charge density at this height). Since the lengths of the sums depend on the prefactor $a$ of the Gaussian, the STM picture generation becomes more expensive if the exponent is raised.

In the constant current mode, the same numerical integration is done, but now the z-position of the current point is scanned in addition until the obtained charge density is equal or larger than the proposed value in -dens. In the first grid point (x=y=0), the imaginary tip is started at z=0.8 times the z-length of the unit cell in order to avoid problems with atoms that might be located at the bottom of the cell at z=0. It is thus assumed that the surface slab is located in the lower half of the unit cell! The smeared local charge density is then calculated by the noted formula and if the charge is still slower than the value given by -dens, the calculation is repeated with the tip lowered by a third of the grid point distance along z in the PARCHG file. This is done until the charge density is equal or larger than the given value or the z-coordinate has reached zero. The corresponding z-value is noted as STM height in the current grid point. Then, the tip is moved to the next x-value. In order to save calculation time, the tip is now not started at z=0.8 times unit cell length, but instead increased only by one Angstroms, assuming that the isosurface of the charge density is sufficiently smooth. Then, the tip is lowered again in small steps until the charge density is high enough, the final z-value is stored. This procedure is repeated until the whole x-y (or more general a-b) grid has been scanned and the STM images has been obtained.

Core level (shift) calculations

This option manages core level shift (CLS) calculations with VASP of systems containing many atoms that shall be evaluated at once. Initial state core level energies can be calculated straightforwardly for all orbitals of all atoms in a system at once, but final state energies need to be calculated for each atom and orbital at once, since the overall electron density needs to be relaxed. This requires a lot of manual work. This program seeks to automatize everything, such that you can calculate initial state and final state energies as well as final state effects of many atoms (all of a certain element) in a system at once.

Setup mode

The setup mode prepares the folders with the global initial state calculation and one final state calculation per selected atom. In the current folder, all input files for a VASP final state calculation need to be present:

• POSCAR (without separation of one atom to the end, as usual for FS)
• POTCAR
• KPOINTS
• INCAR: In this file, all commands for a CLS calculation need to be noted (see section).

The orbital that shall be investigated in each atom of the chosen element needs to be noted there, for example:

ICORELEVEL = 2
CLNT = 4
CLN = 4
CLL = 3
CLZ = 0.5

Here, the 4f electron of the fourth species shall be excited. Note that the species number (CLNT) needs to be one larger than the number of species (elements) in the given POSCAR, since the excited atom in each FS run will be defined as a new separate species. So in this example, the POSCAR and the POTCAR files contain three different elements.

Further, you can prepare a slurm_script file to be used for the IS and each FS calculation.

Then, invoke the program, by specifying the atoms of which element in the system shall be exicted, for example Platinum:

analyze_dft -cls -setup -element=Pt

Now, a folder IS containing the initial state input as well as a number of folders with numbers (e.g., 65,66,67,68) are created. In the enumerated folders, the FS calculations of each atom of the chosen element are listed. In the example, atoms 65, 66, 67, 68 would be all Pt atoms in the system. Now, start all calculations manually.

Else, you can invoke the -slurm option to start all calculations automatically, if you have slurm on your system and prepared a slurm_script in the main folder.

analyze_dft -cls -setup -element=Pt -slurm

Evaluation mode

After you prepared and ran the CLS calculations with the setup mode, simply execute the following command in the main folder:

analyze_dft -cls -eval

Now, the initial state and all final state calculations are evaluated at once. The core level energies of the chosen orbital of the chosen elements are read in, corrected by the Fermi level and further referenced to the pure elemental state by comparing the energy to a predefined list of standard values such that they can be compared directly to experimental values.

Then, the FS, IS and FS-IS energies are plotted as pseudo-spectra, where all energies are broadened by a Gaussian. These spectra are stored in the files plot_fs.dat, plot_is.dat and plot_fs-is.dat.

Finally, the different core level shifts can be visualized spatially with the files show_fs.pdb, show_is.pdb and show_fs-is.pdb. Further, the ranges of energies for each type (FS, IS and FS-IS) are printed, for example:

- show_fs.pdb : Final state values (range:   69.45026  to   70.19829)
- show_is.pdb : Initial state values (range:   69.37304  to   69.86604)
- show_fs-is.pdb : Final state effect values (range:    0.07722  to    0.33225)

Open one of the files with VMD and go to Graphics -> Representations -> Coloring Method -> Occupancy. Now, the atoms will be colored according to for example their FS values.

If you have more than one element in your system, all atoms belonging to a different element were assigned values of zero. Therefore, the range is large, here from 0 to 70.19829, leading to all atoms of the chosen elements being colored blue. To change this, generate a new representation only for the chosen element (for example: name Pt in the representation specifier) and correct the range of energies by going to Trajectory in the Representations window and change there the values in the Color Scale Data Range to those printed out by the script (for example 69.45026 to 70.19829 for the FS pdb file).