Plotting charge density in Ba3Ta2ZnO9

Table of Contents

• Command Summary
• Preliminaries
• Tutorial
  • Plotting the effective potential

This tutorial shows how to extract the charge density from Questaal, in Ba₃Ta₂ZnO₉. It uses a legacy input file that was constructed by hand.

Command Summary

The following commands are sufficient to run the tutorial.

touch ctrl.bzt.
rm *.bzt
cp ~/lmf/b/./fp/test/ctrl.bzt ~/lmf/b/./fp/test/site.bzt .
lmfa bzt
mpirun -n 9 lmf bzt -vnit=30
lmf bzt --rs=0 -vnit=0 --wden:atpos
lmf bzt --rs=0 -vnit=0 --wden:core=2:l1=1,0,0,37:l2=0,1,0,37:fn=atrho
lmf bzt --wden:core=2:l1=1,0,0,37:l2=0,1,0,37

Preliminaries

This tutorial assumes you have cloned and built the lmf repository (located here). For the purpose of demonstration, ~/lmf will refer to the location of the cloned repository (source directory). In practice, this directory can be named differently. It is assumed you have moved the binaries to a directory in your path.

This tutorial runs through a test case taken from the Questaal repository, which you can run as

~/lmf/b/./fp/test/test.fp --quiet --poszer --MPIK bzt

Here we explain and interpret the steps.

This test uses a legacy input file that was constructed by hand. If you want to create a newer looking one, use the site2init utility to create an init file,

site2init site.dat > init.dat

and follow standard procedures to make a ctrl file. Be advised that the automatic procedure will make a somewhat different basis set and it will not generate tags needed for this tutorial, e.g. HAM_GMAX needs to be replaced with HAM_FTMESH.

The tutorial assumes you will run in MPI mode, and uses prefix mpirun -n 9 for MPI jobs. To run in serial mode, just skip the prefix.

At present, there is no capability to interpolate the smoothed density to an arbitrary plane, so you are restricted to choosing a plane that passes through uniform mesh points. (The interstitial density is stored on a uniform grid). To make this translation a little easier, you can run lmf with switch --wden~atpos. This doesn’t generate a density but prints out a table indicating where the nuclei are as multiples of the mesh points.

Tutorial

Start in a fresh working directory and copy the input files needed to run this calculation from the distribution:

 cp ~/lmf/fp/test/{ctrl.bzt,site.bzt} .

Make the free atom density and the self-consistent density

 lmfa bzt
 mpirun -n 9 lmf bzt -vnit=30

Get information about the site positions as multiples of the uniform mesh

lmf bzt --rs=0 -vnit=0 --wden:atpos

You should see a table like this. Note the 36×36×44 mesh of points was specified in the EXPRESS category.

 Site positions as (possibly fractional) indices of 36x36x44 smooth mesh
 Site  p1    p2    p3      Site  p1    p2    p3      Site  p1    p2    p3
   1   0.0   0.0   0.0  |    6   0.0   0.0  22.0  |   11  29.8  23.7  14.3
   2  12.0  24.0  14.9  |    7   0.0  18.0   0.0  |   12  29.8   6.2  14.3
   3  24.0  12.0  29.1  |    8  18.0  18.0   0.0  |   13  12.3   6.2  14.3
   4  12.0  24.0  36.3  |    9  18.0   0.0   0.0  |   14  23.7  29.8  29.7
   5  24.0  12.0   7.7  |   10   6.2  12.3  29.7  |   15   6.2  29.8  29.7

Note the lattice vectors are hexagonal, as can be seen from out put of, e.g., lmchk

                Plat                                  Qlat
   0.000000  -1.000000   0.000000        0.577350  -1.000000   0.000000
   0.866025   0.500000   0.000000        1.154701   0.000000   0.000000
   0.000000   0.000000   1.227430        0.000000   0.000000   0.814710
   alat = 10.926401  Cell vol = 1386.624179

The following writes writes out the density generated by Mattheis construction (superposition of free atomic densities) in the basal plane and passing through the origin, to file atrho.bzt.

$ lmf bzt --rs=0 -vnit=0 --wden:core=2:l1=1,0,0,37:l2=0,1,0,37:fn=atrho

The syntax of --wden is described here; in this context they mean the following:

• core=2    tells lmf to exclude core densities.
• ~l1=1,0,0,37  specifies the first axis for the 2D mesh (1,0,0) corresponds to (0,-1,0) in Cartesian coordinates (see Plat above).
  It selects 37 points because the mesh along P1 consists of 36 points. Thus points 1 and 37 are equivalent.
• ~l1=0,1,0,37  specifies the second axis for the 2D mesh, which is $(\sqrt{3}/2,1/2,0)$. 37 points are specified because P2 also has 36 points.
• fn=atrho   specifies the file name for file I/O.

When writing file atrho.bzt, lmf prints the following:

 ioden : local densities + envelope density, Qloc=19.189235  Q=111.999997
 Writing smooth density to file atrho : origin at (0,0,0).
 Origin, in cartesian coordinates  0.000000 0.000000 0.000000
 v1: (36+1 pts) = (1,0,0)(p1,p2,p3) = (0.000000,-1.000000,0.000000) l=1.000000
 v2: (36+1 pts) = (0,1,0)(p1,p2,p3) = (0.866025,0.500000,0.000000) l=1.000000
 Angle between vectors : 2.094395 radians = 120 deg

 Unit normal vector (0.000000,0.000000,1.000000).  Origin lines in plane 0
 NN planes connected by (0,0,1)(p1,p2,p3) = (0.000000,0.000000,0.027896)
 Planes separated by 0.0278961, 44 planes to shortest period, 0 planes in normal direction

v1 and v2 are the axes for the plane: they are displayed as multiples of Plat and also in Cartesian coordinates. Note the angle between v1 and v2 is 120^o.

Next we write the self-consistent density in the same plane, this time to the default file nane smrho.bzt.

lmf bzt --wden:core=2:l1=1,0,0,37:l2=0,1,0,37

Finally let’s draw a contour plot of the charge density. Rather than plot the density itself, we can plot the shift in the density relative to the Mattheis construction: this shows where the bonding goes.

mcx smrho.bzt atrho.bzt -- > rho_diff.bzt

To draw a figure, copy the following box to file plot.rho

fplot 

# ... for (sqrt(3)/2,1/2,0), (0,-1,0)
 -frme:theta=2*pi/3 .4,1.3,0,.9

  -lt 1,bold=3,col=.00,.1,1.,fill=0 -con  .001 -qr rho_diff.bzt
  -lt 1,bold=3,col=.25,.1,.7,fill=0 -con  .002 -qr rho_diff.bzt
  -lt 1,bold=3,col=.50,.1,.4,fill=0 -con  .005 -qr rho_diff.bzt
  -lt 1,bold=3,col=.75,.1,.1,fill=0 -con  .010 -qr rho_diff.bzt

  -lt 3,bold=5,col=.1,.00,1.0,fill=0 -con -.001 -qr rho_diff.bzt
  -lt 3,bold=5,col=.1,.25,.75,fill=0 -con -.002 -qr rho_diff.bzt
  -lt 3,bold=5,col=.1,.50,.50,fill=0 -con -.005 -qr rho_diff.bzt
  -lt 3,bold=5,col=.1,.75,.25,fill=0 -con -.010 -qr rho_diff.bzt
  -lt 3,bold=5,col=.1,1.0,.00,fill=0 -con -.020 -qr rho_diff.bzt

Then run

fplot -f plot.rho

fplot generates file fplot.ps. Use your favorite postscript viewer to render a drawing of fplot.ps, e.g. gs fplot.ps .

Solid contours correspond to positive densities (actually density differences), varying from blue for most intense, to red for least intense. while dashed contours are for negative densities, varying from blue to green. It is apparent that charge is piling up in a rim around the O atom, and being taken away from a rim around the Ba, consistent with charge transfer from cation to anion in a polar compound. Click below to view the figure fplot should generate.

In units of the lattice vectors, the top right corner is (1,1,0). Atom 8 (an O atom) is located at (1/2,1/2,0) in these units, which lies at the centre of the figure. A Ba sits at the bottom left corner. You can readily see this by invoking

lmchk bzt --pr50

You should see

 site spec        pos (Cartesian coordinates)         pos (multiples of plat)
   1  Ba2       0.000000   0.000000   0.000000    0.000000   0.000000   0.000000
   ...
   8  O1       -0.433013   0.250000   0.000000   -0.500000  -0.500000   0.000000

Plotting the effective potential

You can, instead of plotting the density, plot the (pseudo) effective potential. The true one has singularities at each nucleus; to smooth it out, an effective one is written.

lmf bzt --rs=0 --wden:core=2:l1=1,0,0,37:l2=0,1,0,37:fn=atrho:veff
lmf bzt --wden:core=2:l1=1,0,0,37:l2=0,1,0,37:veff
mcx smrho.bzt atrho.bzt -- > v_diff.bzt

To draw a picture, cut and paste the box below to file plot.pot

fplot

# ... for (sqrt(3)/2,1/2,0), (0,-1,0)
 -frme:theta=2*pi/3 .4,1.3,0,.9

  -lt 1,bold=3,col=.00,.1,1.,fill=0 -con  .01 -qr v_diff.bzt
  -lt 1,bold=3,col=.20,.1,1.,fill=0 -con  .02 -qr v_diff.bzt
  -lt 1,bold=3,col=.40,.1,1.,fill=0 -con  .05 -qr v_diff.bzt
  -lt 1,bold=3,col=.60,.1,1.,fill=0 -con  .1 -qr v_diff.bzt
  -lt 1,bold=3,col=.80,.1,.7,fill=0 -con  .2 -qr v_diff.bzt
  -lt 1,bold=3,col=.99,.1,.4,fill=0 -con  .5 -qr v_diff.bzt

  -lt 3,bold=5,col=.1,.00,1.0,fill=0 -con -.01 -qr v_diff.bzt
  -lt 3,bold=5,col=.1,.25,1.0,fill=0 -con -.02 -qr v_diff.bzt
  -lt 3,bold=5,col=.1,.50,1.0,fill=0 -con -.05 -qr v_diff.bzt

Then run

fplot -f plot.pot

Use your favorite postscript viewer to render a drawing of fplot.ps, e.g. gs fplot.ps.

Questions or Comments