Questaal Home

Plotting charge densities

Table of Contents

Command Summary

The following commands are sufficient to run the tutorial.


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 test case demonstrations generation of charge density in a plane:

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

This tutorial goes through that demonstration, with some additional explanation.

It assumes you will run in MPI mode, and uses prefix mpirun -n 9 for MPI jobs. To run in scalar mode, you can 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.


Start in a fresh working directory of your choice (better not use use the directory where you compiled the code).

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:

 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 . 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 120o.

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


# ... 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
$ open

(Use your favorite postscript viewer for open).

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.

Example 1

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 -- > rho_diff.bzt

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


# ... 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 rho_diff.bzt
  -lt 1,bold=3,col=.20,.1,1.,fill=0 -con  .02 -qr rho_diff.bzt
  -lt 1,bold=3,col=.40,.1,1.,fill=0 -con  .05 -qr rho_diff.bzt
  -lt 1,bold=3,col=.60,.1,1.,fill=0 -con  .1 -qr rho_diff.bzt
  -lt 1,bold=3,col=.80,.1,.7,fill=0 -con  .2 -qr rho_diff.bzt
  -lt 1,bold=3,col=.99,.1,.4,fill=0 -con  .5 -qr rho_diff.bzt

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

Then run

$ fplot -f plot.pot
$ open

Example 2

Questions or Comments

This thread is intended for questions or comments about the tutorial page 'Plotting charge densities' located here.

Feel free to ask questions (or answer!) about the physics found on the page, to get help with the Questaal commands listed or for guidance on how to use Questaal for a similar application that the page did not cover.

For questions about physics or Questaal usage not related to this page, have a look around the site as you may find someone has already asked your question or you can post it yourself in the general questions category.

If you spot a bug, an error or have a suggestion for improvement, please use the GitLab issues page instead.

Anyone is welcome to contribute!