Post-processing CONQUEST output

Introduction

The utility PostProcessCQ allows users to post-process the output of a CONQUEST calculation, to produce structure files, densities of states, and charge density, band densities and STM images as CUBE files (which can be read by the freely available VESTA code).

There are a number of different analyses which can be performed: coordinate conversion (to formats which can be plotted); conversion of total charge density to CUBE file format; production of band-resolved (optionally k-point resolved) densities in CUBE file format; simple Tersoff-Hamann STM simulation; and calculation of densities of states, including projected DOS. You should ensure that all the files produced during the CONQUEST run are available for the post-processing (including eigenvalues.dat, chden.NNN, make_blk.dat or hilbert_make_blk.dat and ProcessNNNNNNNWF.dat and ProcessSijNNNNNNNWF.dat as applicable) as well as the input files.

Note that the utility reads the Conquest_input file, taking some flags from the CONQUEST run that generated the output, and some utility-specific flags that are detailed below.

Projected DOS, band density and STM simulation are available for calculations using primitive PAOs or multi-site support functions (MSSFs), but not for calculations using blip functions.

Go to top.

Coordinate conversion

Set Process.Job coo to output a coordinate file for further processing or plotting. The utility will read the file specified by Process.Coordinates (which defaults to the file specified by IO.Coordinates). The output format is selected by specifying the Process.CoordFormat tag. The default output format is XYZ (which adds a .xyz suffix to the file name) using xyz. The CASTEP .cell output format can also be selected using cell. We plan to expand this conversion to other formats in the future.

Note that for a structural relaxation or molecular dynamics calculation, if you do not specify Process.Coordinates then the IO.Coordinates file, which will be converted, will be the input structure, not the output structure. Parameters that can be set are:

Process.Coordinates string (default: IO.Coordinates value)
Process.CoordFormat string (default: xyz; options: xyz, cell)

Go to top.

Charge density

Setting Process.Job to cha, chg or den will convert the files chden.NNN which are written by CONQUEST to a cube file. The processing will use the files chden.NNN, Conquest_input and hilbert_make_blk.dat or raster_make_blk.dat. Parameters that can be set include:

Process.ChargeStub string (default: chden)

The ChargeStub simply defines the filename which will be read, and used for output.

Note that to output the chden.NNN files from CONQUEST, you must set the flag IO.DumpChargeDensity T in the CONQUEST run.

Go to top.

Band density

Setting Process.Job to ban produces band densities from wave function coefficients output by CONQUEST. The non-self-consistent, dense k-mesh run must read in a well-converged charge density and must have the following tag set:

IO.outputWF T

By default, all bands will be produced. Restricting the set of bands whose coefficients are output can be done with either an energy range:

IO.WFRangeRelative T/F
IO.min_wf_E real (Ha)
IO.max_wf_E real (Ha)

or with a list of bands:

IO.maxnoWF n

%block WaveFunctionsOut
n entries, each a band number
%endblock

The wavefunction range can be relative to the Fermi level (IO.WFRangeRelative T) otherwise it is absolute. Either of these will produce a file containing all eigenvalues at all k-points (eigenvalues.dat) and a series of files containing the wavefunction expansion coefficients for the selected bands (ProcessNNNNNNNWF.dat). These files are output as binary (unformatted) by default (this can be changed by setting IO.MatrixFile.BinaryFormat F before the CONQUEST run) and will be read using the same format (it is important to check this!).

From these wavefunction coefficient files, band densities can be produced in post-processing, using similar tags; either a range:

Process.min_wf_E real (Ha)
Process.max_wf_E real (Ha)
Process.WFRangeRelative T/F

or an explicit list of bands:

Process.noWF n

%block WaveFunctionsProcess
n entries, each a band number
%endblock

Note that the bands to be processed must be a subset of the bands output by CONQUEST. The bands can be output summed over k-points, or at individual k-points, by setting Process.outputWF_by_kpoint to F or T respectively.

The output log from post-processing will contain lines of the form "Integral of band n with energy E, is <float>" (Process.outputWF_by_kpoint T) or "Integral of band n at kp N is <float>". A sufficiently converged integral will be around 1.00 for every k-point and band.

Go to top.

Tersoff-Hamann STM simulation

Setting Process.Job ter will use a very simple Tersoff-Hamann approach to STM simulation, summing over band densities between the Fermi level and the bias voltage (this is often surprisingly accurate). The following parameters can be set:

STM.BiasVoltage    real (eV)
STM.FermiOffset    real (eV)
Process.MinZ       real (Bohr)
Process.MaxZ       real (Bohr)
Process.RootFile   string (default: STM)

The FermiOffset tag allows the user to shift the Fermi level (to simulate charging or an external field). The height of the simulation cell in which the STM image is calculated is set by the MinZ and MaxZ tags, and the filename by the RootFile tag.

Go to top.

Density of states (DOS)

Setting Process.Job dos will produce a total density of states (DOS) for the system, using the eigenvalues.dat file output by a CONQUEST diagonalisation calculation. The following parameters can be set:

Process.min_DOS_E       real    (Ha, default lowest eigenvalue)
Process.max_DOS_E       real    (Ha, default highest eigenvalue)
Process.WFRangeRelative T/F
Process.sigma_DOS       real    (Ha, default 0.001)
Process.n_DOS           integer (default 1001)

The limits for the DOS are set by the first two parameters (note that CONQUEST will output all eigenvalues). These limits can be given as absolute energies, or relative to the Fermi level, depending on whether WFRangeRelative is set to F or T respectively. The broadening applied to each state is set by sigma_DOS, while the number of bins is set by n_DOS. The integrated DOS is also calculated; the user can choose whether this is the total integrated DOS (i.e. from the lowest eigenvalue, regardless of the lower limit for DOS) or just the local integrated DOS (i.e. over the interval specified for the DOS) by setting Process.TotalIntegratedDOS to T or F, respectively.

If the user does not specify limits on the DOS, then the range will be expanded slightly so that the broadened peaks fit completely within the DOS range. This behaviour can be controlled with an additional parameter Process.ExpandRange T/F (the default behaviour is false if the user sets any limits, otherwise true).

We recommend that, for accurate DOS, CONQUEST should be run non-self-consistently with a very high k-point density, after reading in a well-converged input charge density: set minE.SelfConsistent F and General.LoadRho T (which will require that the converged charge density is written out by CONQUEST by setting IO.DumpChargeDensity T).

Go to top.

Atom-projected DOS

Setting Process.Job pdos will produce a total density of states as above, as well as the density of states projected onto the individual atoms. Given support functions \(\phi_{i\alpha}(\mathbf{r})\) which are the basis functions of the Kohn-Sham eigenstates \(\psi_{n}(\mathbf{r}) = \sum_{i\alpha} c^{n}_{i\alpha}\phi_{i\alpha}(\mathbf{r})\), then the projection of a given state, \(n\), onto an atom \(i\) can be written as \(\sum_{\alpha j\beta} c^{n}_{i\alpha} S_{i\alpha,j\beta}c^{n\mathbf{k}}_{j\beta}\). The projected DOS is constructed using these projections.

For calculations using pseudo-atomic orbitals (PAOs), either directly or through multi-site support functions, the atom-projected DOS can be further resolved by angular momentum (either just \(l\) or both \(l\) and \(m\)). For MSSF calculations, CONQUEST outputs the required coefficients in the underlying PAO basis. Projected-DOS post-processing is not available for calculations using blip functions. To output the necessary coefficients, a static diagonalisation calculation must be performed with:

IO.write_proj_DOS T

As for the DOS, very high Brillouin zone sampling is required for accurate projected DOS, which is most efficiently generated using a converged charge density and a non-self-consistent calculation with much higher k-point density. CONQUEST will produce the wavefunction files (ProcessNNNNNNNWF.dat and ProcessSijNNNNNNNWF.dat) as binary (unformatted) by default (change using the flag IO.MatrixFile.BinaryFormat F).

Once the files have been generated by CONQUEST, the output can be processed by setting the output tag:

Process.Job pdos

This is all that is needed for the simplest output. The number of bins and smearing of the peaks can be set using:

Process.sigma_DOS 0.002
Process.n_DOS 10001

To resolve the DOS by angular momentum as well as by atom, then the following flags can be set:

Process.pDOS_l_resolved T
Process.pDOS_lm_resolved T

Note that only one of these is needed, depending on what level of resolution is required. Angular-momentum resolution is available for both primitive-PAO and MSSF calculations, but not for calculations using blip functions.

The energy range for the projected DOS can also be specified:

Process.min_DOS_E -0.35
Process.max_DOS_E  0.35
Process.WFRangeRelative T

where the final tag specifies that the limits on the DOS are given relative to the Fermi level.

If you only want to produce pDOS for a few atoms, then you can set the variable Process.n_atoms_pDOS and list the atoms you want in the block pDOS_atoms:

Process.n_atoms_pDOS 2
%block pDOS_atoms
1
12
%endblock

Go to top.

Rotated pDOS

CONQUEST, by default, projects the orbitals along the simulation cell axes. Sometimes it is convenient to be able to project onto a coordinate system defined by an atom’s local environment, e.g. bonds, to facilitate chemical analysis. CONQUEST supports rotating the wavefunction coefficients either by inputting the basis of the final coordinate system, using the algorithms described by Maintz et al. and Romanowski et al. [PP1, PP2], or by supplying active Euler angles in the extrinsic \(zyz\) convention.

There are four operating modes:

  1. Supply a set of axes to apply to a set of (or all) atoms (mode 0)

  2. Supply 3 Euler angles to apply to a set of (or all) atoms (mode 1)

  3. Supply specific atoms and information about their neighbours to construct local axes for each atom (mode 2)

  4. Supply specific atoms each with a rotation axis and angle (mode 3)

To enable any form of pDOS rotation, the flag Process.RotatePDOS T must be set in Conquest_input. It defaults to F, which will use the standard projection onto the simulation cell axes.

To use mode 0, have the following snippet in Conquest_input:

Process.RotatePDOS T (default: F)
Process.RotatePDOSMode 0 (default)
%block pDOSAxes
x1 y1 z1 ! new x
x2 y2 z2 ! new y
x3 y3 z3 ! new z
%endblock
Process.n_atoms_pDOS 2
%block pDOS_atoms
1
12
%endblock

The above snippet says to use the user-input axes defined in the pDOSAxes block. The first line defines the new \(x\) direction along vector \((x_1, y_1, z_1)\), then \(y\) along \((x_2, y_2, z_2)\) and \(z\) along \((x_3, y_3, z_3)\). In this mode, all atoms will have this rotation applied to them by default. The pDOS_atoms block states which of these are written out to disk. The code will abort if the vectors provided do not form a right-handed coordinate system, or if the vectors are not orthogonal. It is recommended to supply the directions as vectors of integers - the code will take care of normalisation.

If Euler angles are desired, set Process.RotatePDOSMode 1 and in Conquest_input use the following lines:

Process.RotatePDOS T
Process.RotatePDOSMode 1
Process.RotatePDOS.NumAtoms int (required)
Process.RotatePDOSAngle deg (default: deg, optional: rad)
%block pDOSEuler
A alpha beta gamma !n entries, equal to Process.RotatePDOS.NumAtoms
%endblock pDOSEuler

Here, A is the atom number to assign the Euler angles to. The default angle unit is degrees (Process.RotatePDOSAngle deg), and units cannot be mixed. In the extrinsic \(zyz\) convention, three consecutive rotations happen about a fixed set of coordinates, i.e. the CONQUEST simulation-cell axes (CONQUEST only supports orthorhombic cells): first about the \(z\) axis by \(\gamma\), then the \(y\) axis by \(\beta\), and then the \(z\) axis by \(\alpha\), in a right-handed sense (looking down any of the axes towards the origin, the positive direction of rotation is anticlockwise).

Determining local axes is possible depending on local geometries.

Process.RotatePDOS T
Process.RotatePDOSMode 2
Process.RotatePDOS.NumAtoms int (required)
%block pDOSNeighbours
A B C D  !n entries, same number as in Process.RotatePDOS.NumAtoms
%endblock pDOSNeighbours

where A, B, C, D are integers with different conditions:

  • A: the atom number to perform the rotation for

  • B: if it is 0, assume square planar geometry; if 1, assume octahedral. This controls nearest-neighbour searching for the nearest four or six neighbours, respectively.

  • C: if -1, set the principal direction along the shortest bond. If 0, set along the longest bond. If it is a positive integer, it is interpreted as an atom number which must be a neighbour of the atom specified in A. All runs will output the atom neighbours, so setting this to 0 or -1 as a first run is recommended unless the neighbours are known beforehand

  • D: if 0, the second direction is chosen by the bond which changes the least under projection to the plane defined by the bond vector calculated from the neighbour specified in C. If D is a positive integer, then it is interpreted as an atom number which must be a neighbour of the atom specified in A. A warning will be displayed if the direction towards neighbour D is determined to deviate significantly away from orthogonal.

If the rotation axis and angle are known beforehand, using mode 3 can be convenient.

Process.RotatePDOS T
Process.RotatePDOSMode 3
Process.RotatePDOS.NumAtoms int (required)
%block pDOSAxisAngle
atom_number x y z angle  !n entries, same number as in Process.RotatePDOS.NumAtoms
%endblock

The atom number that this defines a local axis for must be specified. Then it is followed by three floats specifying the rotation axis in the basis of the simulation cell coordinate system. Finally, this is followed by a rotation angle which is positive when looking down the rotation axis and rotating anticlockwise.

For modes 1, 2 and 3, setting the block pDOS_atoms is optional since only the specified atoms will be rotated whilst the rest will be in the simulation cell axes.

For users interested in verbose output, set Process.RotatePDOSDebug T (default F) in Conquest_input which will output detailed information about rotation matrices and orbital weights.

Go to top.

Band structure

The band structure of a material can be generated by CONQUEST by performing a non-self-consistent calculation, after reading a well-converged charge density: set minE.SelfConsistent F and General.LoadRho T (remember that to write a converged charge density from CONQUEST you set IO.DumpChargeDensity T). The k-points required can be specified as lines of points in k-space. With Diag.MPMesh F, setting Diag.KspaceLines T enables this, while the number of lines (e.g. Gamma to L; L to X; would be two lines) is set with Diag.NumKptLines and the number of points along a line with Diag.NumKpts. The k-point lines themselves are set with a block labelled Diag.KpointLines which should have two entries (starting and finishing k-points) for each k-point line. (In constructing the k-point list, CONQUEST will automatically remove any duplicate points, so that the output can be plotted smoothly.) So to create a bandstructure from X-\(\Gamma\)-L-X (3 lines: X-\(\Gamma\); \(\Gamma\)-L; L-X) with 11 points in each line, you would use the following input:

Diag.KspaceLines T
Diag.NumKptLines 3
Diag.NumKpts 11
%block Diag.KpointLines
0.5 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.5 0.5 0.5
0.5 0.5 0.5
0.5 0.0 0.0
%endblock

After running CONQUEST, setting Process.Job bst and running the post-processing will read the resulting eigenvalues.dat file, and produce a file BandStructure.dat. The x-axis will be the k-point index by default, but specifying Process.BandStrucAxis (taking value n for index, x, y or z for a single direction in k-space, or a to give all k-point coordinates) will allow you to control this. Limits on the energies to select the bands produced can be set with Process.min_DOS_E and Process.max_DOS_E.

Go to top.

[PP1]

S. Maintz, M. Esser, and R. Dronskowski. Efficient rotation of local basis functions using real spherical harmonics. Acta Physica Polonica B, 47(4):1165, 1 2016. URL: https://doi.org/10.5506/aphyspolb.47.1165, doi:10.5506/aphyspolb.47.1165.

[PP2]

Zbigniew Romanowski, Stanisław Krukowski, and Abraham F. Jalbout. Transformation of Real Spherical Harmonics under Rotations. Acta Physica Polonica B, 39(8):1985–1991, 8 2008. URL: https://www.actaphys.uj.edu.pl/fulltext?series=Reg&vol=39&page=1985.