Behaviour of classes Path (lattice.py) and YambopyBandStructure (bandstructure.py)

[palful]

It is not clear if these classes actually produce correct band structures for non-hexagonal systems. Both in finding correctly the path and in calculating the relative distances between kpoints. In general, the Path class is unreliable for k-paths in non-hexagonal systems and there is too much redundancy in the code, with the function expand_kpoints appearing in three different places (lattice.py, YamboLatticeDB, YamboSaveDB).

This must be corrected and streamline. I report various issues below, some taken from previous issues.

[palful]

One issue in particular is that it seems to experience problems in finding the correct relative distances between the points of the given paths in a rectangular case: the k-point of the X and Y points are given in crystal coordinates as [0.5,0.0,0.0] and [0.0,0.5,0.0], while the lattice vectors a and b have different lenghts. However, the segments GX and GW displayed in the band structure have the same length.

[palful]

For some paths in the BZ, get_path will fail to work properly. In particular, with Marco we verified that the GXLW path in fcc systems is not correctly generated.

Oddly enough, this only applies if the band structure is NOT interpolated. Alejandro verified that interpolation works fine.

Accurate testing of get_path in both QE (data-file_schema) and yambo (SAVE/ns.db1) is in order.

Proposal: automatic testing by travis with all the paths in https://en.wikipedia.org/wiki/Brillouin_zone?

[palful]

Pedro proposed to add some default hardcoded paths in the BZ when plotting band structures: in this way they don't need to be given explicitly by the user every time.

List of paths for every crystal system here: https://en.wikipedia.org/wiki/Brillouin_zone

[palful]

For calculations in rombohedral lattices yambopy does not get enough points along the path. It works if running the interpolation mode.

[jmcastelo]

For some paths in the BZ, get_path will fail to work properly. In particular, with Marco we verified that the GXLW path in fcc systems is not correctly generated.

The new BrillouinZone class (perhaps to be renamed Path and substitute it) constructs paths in the BZ according to crystal system. Check it in my branch! When constructing the GXLW path with it for a FCC crystal and feeding it to YamboLatticeDB's get_path, the 3D plot of the expanded k-points, the path points and those co-linear with the path shows that (apparently) no expanded point which is co-linear is left out.

I had to redefine the FCC lattice vectors according to QE's criterion for it to work. When using Curtarolo's, the points selected by get_path did not coincide with the path, but were located in an equivalent area of the BZ. If someone can throw light on this, it will be welcome.

[palful]

Following a discussion with @corentinmorice1, I can pinpoint a couple of issues of the old Path which should be avoided in the new implementation. It also ties into José's comment above.

Issue 1: when we do a band plot, we typically want to plot the shortest distances between equivalent high-symmetry points. Let us use the example of the FCC lattice. Suppose that we are doing the path GXLW: we want all these points to be specified in the same quadrant of the BZ, despite the fact that you have many possible equivalent coordinates. For example, you have 8 equivalent L points (at the center of each hexagonal face), but we want to specify the coordinates of the one closest to the other points, not the one on the opposite face. Otherwise, the XL and LW distances will be much longer. The actual location of a specific L-point, say, the one with coordinates (1/2,1/2,1/2), as opposed to the one with (-1/2,-1/2,-1/2), will also depend on the orientation of the lattice vectors.
Action: it could be desirable that for plotting band structures, the user could specify any one of the equivalent high-symmetry points, but yambopy will always compute distances between the closest equivalent points.

Issue 2: not all the equivalent points are sampled in the kpoint meshes, only half of them. So for example, if a certain L point is sampled, the equivalent one on the opposite side (i.e. -L which is still L) will be absent from the sampling. (This is because it is reachable from L by summing an integer-coordinate reciprocal lattice vector (-1,-1,-1), therefore including both -L and L would be a double counting). At present, sometimes Path fails when the coordinates of a non-sampled high-symmetry points are specified by the user instead of the coordinates of the sampled one.
Action: we could make sure that whether or not the high-symmetry point is sampled, yambopy always calculates the correct paths.

[palful]

Addendum following the discussion with Corentin, in case it may be relevant for the new patch. In an hexagonal system the Path distances from the middle to the edge of the BZ (Gamma-A, M-L, K-H) come out with the wrong length in the current version, as can be seen from this plot comparing yambopy (right) to ASE (left) from the same calculation. In this particular example, there are no intermediate points in the Path along the k_z direction, only the high-symmetry ones. Maybe this is the cause of the problem.

[jmcastelo]

Issue 1 above by @palful makes me think of the function expand_kpoints of YamboSaveDB and YamboLatticeDB classes. This function: "Takes a list of k-points and symmetry operations and returns the full Brillouin zone with the corresponding index in the irreducible Brillouin zone." I haven't thought much about how to implement the action associated with this issue, but perhaps this function will be helpful. Maybe expanding the high-symmetry points?

As for Issue 2, I don't understand what is meant by "sampling". Is it the mesh that QE computes?

[jmcastelo]

There's another issue with the new Path class I am developing. Curtarolo's criterion to define the primitive lattice vectors of the crystal types sometimes differs from QE's. This implies that the reciprocal lattice vectors may be different and so the high-symmetry points. If the user makes a QE calculation specifying an Bravais lattice with ibrav (different from 0), currently there's a discrepancy between the high-symmetry points of both criteria. Since the new Path uses the former criterion to define paths in the BZ, the result is that get_path (now located outside its original classes) doesn't find collinear points of the expanded k-mesh corresponding to the QE calculation.

Maybe there's a way to transform from a representation of the high-symmetry points on one criterion to the other. The user could specify that the desired path corresponds to QE's criterion and the class would get the right points from Curtarolo's.

I don't know how to define high-symmetry points with the reciprocal lattice bases that QE uses and associates with ibrav different from 0. So far, I am sticking to ibrav=0 and entering the lattice vectors as Curtarolo defines in his paper.

[palful]

My short comment is that I think that, if we have to adopt a certain choice for lattice vectors, we should use the one of QE since Yambo is run on top of QE and therefore in virtually 100% of the cases we will have that convention.

It may be tedious to find the correct coordinates for the high-symmetry points of choice, but it could be done for example using xcrysden to read a sample input file for every ibrav type, and then using its k-path utility (Alejandro knows how to do it). Or, we could have a look at how ASE does the same thing. Importing ASE as a library in order to use its functionalities is certainly possible.

Below a more elaborate comment.

Here we have an ambiguity since there are many equivalent high-symmetry points. For example, in a 2D hexagonal lattice, we have 6 equivalent M points (in the middle of the hexagon sides) and 6 equivalent K points (at the vertices). Moreover, only half of them are actually contained in the mesh that QE computes, since the others are assigned to adjacent Brillouin zones. To compound this, the coordinates of a specific high-symmetry point (like the M-point on the "top" side of the hexagon) depend on the choice of primitive lattice vectors.

To clarify this, below we have two possible choices of primitive lattice vectors for the 2D hexagonal cell (the vector lengths are not at scale and not placed very precisely). The vectors in red correspond to the convention of Quantum Espresso. The orange points are those sampled in a standard calculation using symmetries, the teal points are the result of the expansion done by yambopy and yambo. The indices of each point in the two meshes (orange - original, teal - expanded) are also shown. Notice also how the sides of the BZ are not fully sampled. In order to obtain the points not sampled, one needs to consider adjacent BZs as well. This is something that I think yambopy does, I was discussing it with Ali recently.

Starting from the expanded mesh, in principle any path connecting any high-symmetry point can be drawn (the standard here is Gamma-M-K-Gamma). This circuit can be done in any part of the BZ, since all the points will be equivalent. For convenience, typically one chooses the path in the same quadrant of the irreducible BZ (orange area), so our Gamma-M-K-Gamma circuit will be on the top right portion of the hexagon most likely.

Below a plot to show the position of some of the equivalent M and K points (some sampled in the first BZ, some not), along with the fact that their coordinates may change depending on the choice of lattice vectors. (The coordinates of the alternative convention, the blue arrows, may not be all correct, I kind of filled them quickly just to illustrate the point).

Clearly, it is important that the user knows the location of the points they are putting in the Path, we can't do that for them. However, if one chooses a default path with default points, we have to make sure that we are doing the correct circuit, so when we are calculationg the GM, MK, KG segments, they should be the shortest possible between all equivalent points.

If we need to use the coordinates relative to a specific choice of primitive lattice vectors, then these should refer to the choice of primitive lattice vectors by QE, which is most likely what the users will have. If they are using ibrav=0, then it's their responsibility to know the location of the points and the Path they want.

[jmcastelo]

Thanks for the detailed and clear comment, @palful !
I'm rewriting the BrillouinZone class to make it QE-compatible.

However, in the classification of BZs given by Curtarolo's paper, the reduced coordinates of some high-symmetry points depend on the lattice parameters, for some crystal lattices. This makes me think that perhaps in those cases it may not be straightforward to define them given just the primitive lattice vectors. This applies, for example, to the body-centred tetragonal lattice.

Moreover, it appears that there are BZ variants (BCT1 and BCT2 for the body-centred tetragonal lattice), whose high-symmetry points reduced coordinates are given by different formulas.

Any idea about how to solve this?

[jmcastelo]

Clearly, it is important that the user knows the location of the points they are putting in the Path, we can't do that for them. However, if one chooses a default path with default points, we have to make sure that we are doing the correct circuit, so when we are calculationg the GM, MK, KG segments, they should be the shortest possible between all equivalent points.

The BrillouinZone class already allows the user to define a path as a set of segments by selecting the labels of the available and predefined high-symmetry points of the crystal lattice of choice.

Example for a BCC: [['Gamma', 'H', 'N', 'Gamma', 'P', 'H'], ['P', 'N']]

Since these points are predefined, the issue with the segments' lengths being the shortest between all equivalent points may be solved by an appropriate selection of their reduced coordinates. One can do as suggested, namely with xcrysden plotting the first BZ corresponding with an ibrav and choosing the high-symmetry points that fulfill the requirement that they are the closest to each other. These reduced coordinates are then hard-coded into BrillouinZone.

Maybe the QE's k-points mesh is in another area of the first BZ than the defined path, but because the mesh is expanded using the symmetries of the crystal, and shifted to adjacent BZs, the function get_path will select the collinear points (if any) without issues (that I can foresee). The resulting collinear points and their indices can be used to select the energies (given now by YamboSaveDB) and passing both path and energies to YambopyBandstructure, the band structure plots can be done. (I modified this last class to accept a BrillouinZone path.)

I am systematically testing all ibravs. In my tests, a default path is chosen for each ibrav, and some QE calculations are performed with meshes of varying k-point numbers (8x8x8, 9x9x9, 12x12x12, etc.) to check if some collinear points are detected. A 3D plot is done to visualize the mesh, path and collinear points, equivalent to the 2D plots above. Finally, a plot of the band structures is done too. If anyone wants me to perform a test on a specific calculation, please send me the input files and I will do it.

[palful]

I have a small update which I addressed with PR #48 .

In order to plot a band structure, one should use cartesian coordinates to preserve the ratio between distances. However, both in QPDB and ExcitonDB reduced coordinates were used. This gave correct results only for the 2D hexagon apparently. The reduced coordinates are used by the Path object (because it's easier for the user to set it with these coordinates) and by the SKW interpolator, which works in rlu. In the code, they were never changed to cc before being given to the BandStructure object. By making this change, the band structures of rectangular lattices are now correct.

I don't know if this was the only issue and there are no problems with get_path after all. But it works at least as a temporary solution. So far I also created unique copies of expand_kpoints and get_path inside yambopy.kpoints, removing them from the various places they previously appeared in.

[palful]

Now this issue and issue #31 (yambopy-ASE interface) are the same.