In theGreatWhiteShark/vasp2R: R package to import and plot results obtained by the DFT program VASP.

About

This vignette shows the basic usage of the vasp2R package. First we will import an example CHG file, do some manipulation and finally plot the result of the VASP calculation.

The VASP files used to create the CHG file supplied by this package are taken from the third Hands-on example of VASP's documentation. I used the first one: the clean Ni surface. Why did I used the CHG and not the CHGCAR? Because it also contains the required information but is smaller in size and I want to keep package's size at a minimum.

Import

First of all you have to be sure of have the vasp2R package installed.

check.validity <- try( library( "vasp2R" ), silent = TRUE )
if ( class( check.validity ) == "try-error" ){
    ## the package is not installed yet
    devtools::install_github( "theGreatWhiteShark/vasp2R" )
    library( vasp2R )
}

During installation the CHG files is copied as one of the package's assets on the hard disk. We will import it into R using the vasp.import function.

data.imported <- vasp.import( system.file( "example/CHG",
                                          package = "vasp2R" ) )

Structure of the imported vasp object

The resulting object data.imported is of class vasp and inherits from the list. It has three named elements: charge, atoms and lattice

charge

The charge element represents the charge density in the three dimensional unit cell. It itself is of class data.frame and contains four named columns: x, y, z, charge.

atoms

The atoms element represents the coordinates and the type of the individual atoms in the unit cell. It is derived from the CONTCAR like entry in the CHGCAR/CHG file. So the type value corresponds to the VASP internal string representation of the periodic table (e.g. "Si" for silicon). It is of class data.frame and consists of four named columns: x, y, z and type

lattice

The lattice element is a data.frame representing the lattice vectors of the unit cell. It contains three named columns: x, y and z and the entries of a specific row correspond to a specific vector.

Manipulation of the vasp object

Reproducing the unit cell

Afterwards we will reproduce the imported VASP results of a single Ni atom with periodic boundary conditions to make it actually look like a surface

data.reproduced <- vasp.reproduce( data.imported, x.rep = seq( -3, 3 ),
                                  y.rep = seq( -3, 3 ) )
## the x.rep, y.rep and z.rep argument of the vasp.reproduce function
## have to be supplied using a numerical vector. In the above example
## the unit cell is going to be reproduced three times in the positive
## and negative direction along both the x and y axis.

Most of the functions in vasp2R can be called with either the whole object of class vasp or with the individual elements. So its also possible to call vasp.reproduce( data.imported$atoms ) to calculate just a reproduced version of the atomic positions.

Rotating the surface

We can also rotate the positions to change our point of view. For such a clean and symmetric surface this is a little bit useless. But for more complex structures it will become quite handy.

data.rotated <- vasp.rotate.cell( data.reproduced, angle = 45* pi/ 180 )

Creating bonds between atoms

Then we establish bonds between the individual atoms (for a better visualization of the results).

data.bonds <- vasp.bonds( data.rotated, distance = 8 )

The vasp.bonds function will create bonds between all atoms closer together than the supplied distance value within the x-y plane.

Further modifications

Another function implemented in vasp2R is vasp.diff. It expects two objects of class vasp as input and creates a new one of the same type containing the charge differences of the first input minus the second one. It also holds all the atoms and lattice entries of the first argument and is therefore a fully fledged vasp object.

Visualization

In principle you can use whatever plotting package you want. But personally I would recommend the usage of the ggplot2 package.

But before we start visualizing anything, how do we plot the charge density? My approach is to slice the three dimensional charge density in individual planes and plot those along with the atoms in the plane.

There are two main problems in visualizing the charge density:

1. The spacing of the charge density values provided by the CHGCAR is not compatible with the one the ggplot2 libraries accept. This actually leads to a huge number of tiny rectangle getting plotted and you won't even be able to see. That's how tiny they gonna be.
2. There are way to much charge density values. If one intends to save the plot as a vector graphic you would end up with several GB even for moderate sized unit cells (simply because they contain hundreds of thousand points)

So instead we create an evenly grid of variable size and calculate the mean charge density in each of those grid boxes.

library( ggplot2 )
library( RColorBrewer ) # For nice colors out of the bxo
color.atom <- "navy"
color.gradient <- "YlOrRd" # a RColorBrewer palette

## Here we will create a plot function taking a z value, create
## a slice at this specific value and plots it along with the atomic
## positions in the x-y plane
vasp.plot <- function( vasp.input, z = 0, grid.point.number = 20 ){
    ## Check out RColorBrewer::brewer.pal.info for different ones

    ## Determine the z grid point closed to the supplied z value
    ## This is not the actual height value but its index in
    ## vasp.input$charge$z!
    plot.height <- which( abs( unique( vasp.input$charge$z ) - z ) == 
                          min( abs( unique( vasp.input$charge$z ) - z ) ) )
    ## Extract just those values of the charge density corresponding to
    ## the determined slice grid point
    slice.charge <-  vasp.input$charge[
        vasp.input$charge[ "z" ] == unique( vasp.input$charge$z )[
                                              plot.height ], ]
    ## Same for the atoms
    plot.atoms <-  vasp.input$atoms[
        vasp.input$atoms[ "z" ] == unique( vasp.input$atoms$z )[
                                             plot.height ], ]

    ## Constructing a grid for plotting the charge density
    ## Boundaries of the charge density
    limits.x <- c( min( slice.charge$x ), max( slice.charge$x ) )
    limits.y <- c( min( slice.charge$y ), max( slice.charge$y ) )

    ## Creating the grid
    plot.charge <- expand.grid( x = seq( limits.x[ 1 ],
                                          limits.x[ 2 ], , 
                                          grid.point.number ),
                                 y = seq( limits.y[ 1 ],
                                         limits.y[ 2 ], ,
                                         grid.point.number ) )

    ## Width of the grid boxes
    grid.width.x <- unique( plot.charge$x )[ 2 ] -
                      unique( plot.charge$x )[ 1 ]
    grid.width.y <- unique( plot.charge$y )[ 2 ] -
                      unique( plot.charge$y )[ 1 ]

    ## Calculate the mean charge density for all points within one grid
    ## box
    plot.charge$charge <- Reduce( c, apply( plot.charge, 1, function( row ){
        charge <- mean( slice.charge$charge[
            abs( slice.charge$x - as.numeric( row[ 1 ] ) ) <=
                grid.width.x &
            abs( slice.charge$y - as.numeric( row[ 2 ] ) ) <=
                grid.width.y ] )
        return( charge )
    } ) )

    ## the actual plotting
    ggplot() +
        ## Plotting the charge density of the slice
        geom_raster( data = plot.charge, interpolate = TRUE,
                  aes( x = x, y = y, fill = charge ) ) +
        ## Plotting the atoms
        geom_point( data = plot.atoms, color = color.atom, stroke = 1.5,
                   aes( x = x, y = y ), shape = 16, size = 4 ) +
        ## Customized the color scale of the charge density
        scale_fill_distiller( palette = color.gradient, direction = 1 ) +
        ## Remove the grid lines in the ggplot plot
        theme_bw() + 
        theme( panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank(), 
              panel.background = element_blank() )

    ## Via this functions additional elements can be introduced to the
    ## plot after calling this function
    return( last_plot() )
}

vasp.plot( data.imported, 0, 20 )

And for the reproduced one

vasp.plot( data.reproduced, 0, 20 )

Pretty good so far. Now lets introduce the bonds to the plotting.

The data.bonds$bonds element contains the point of the beginning and end of all the extracted bonds.

vasp.plot( data.bonds, 0, 20 ) +
    ## extract just the bonds within the considered plane
    geom_segment( data = data.bonds$bonds[
                      data.bonds$bonds$z.begin == 0 &
                      data.bonds$bonds$z.end == 0, ],
                 aes( x = x.begin, y = y.begin, xend = x.end,
                     yend = y.end ), colour = color.atom, size = 1 )

Now we have bonds to the nearest neighbors.

If we want to have additional bonds to the next nearest neighbors we just have to increase the distance argument in the vasp.bonds function. (this time without rotation)

data.plot <- vasp.bonds( data.reproduced, distance = 15 )

vasp.plot( data.plot, 0, 20 ) +
    ## extract just the bonds within the considered plane
    geom_segment( data = data.plot$bonds[
                      data.plot$bonds$z.begin == 0 &
                      data.plot$bonds$z.end == 0, ],
                 aes( x = x.begin, y = y.begin, xend = x.end,
                     yend = y.end ), colour = color.atom, size = 1 )

theGreatWhiteShark/vasp2R documentation built on May 31, 2019, 9:10 a.m.