Plane-wave reader - G-vector ordering and k-dependent basis (AO) size

[kousuke-nakano]

Dear @scemama and @q-posev

Hello. @selim-tunbul and I(@kousuke-nakano ) are implementing a Quantum ESPRESSO (QE) reader for the plane-wave section. We’ve encountered two issues.

Ordering of reciprocal vectors

One challenge is that the ordering of reciprocal-lattice vectors ($G$-vectors) in QE is not guaranteed to match the ordering (or assumptions) in TREXIO. In QE, the SCF output at each $\mathbf{k}$ point is a set of coefficients $c_{i,k,m}$ for each eigenstate $i$ (band) and each $G$-vector index $m$. The $i$-th molecular orbital at $\mathbf{k}$ expands as

$$ \Psi_{i,\mathbf{k}}(\mathbf{r})=\sum_m c_{i,\mathbf{k},m}, e^{,i(\mathbf{k}+\mathbf{G}_m)\cdot \mathbf{r}}. $$

TREXIO can store the $i$ and $k$ indices, but there appears to be no explicit place to store the ordering (or the list) of G-vectors indexed by $m$.

In QE, we can recover the order because the Miller indices associated with each $G$-vector are provided. In TREXIO’s plane-wave section, however, the only relevant parameter seems to be the cutoff energy (e_cut), and no canonical $G$-vector list or ordering is specified.

We see two paths forward:

1. Add an explicit dataset for $G$-vectors (or their order) to TREXIO.
   For example, store the Miller indices $(h,k,l)$ for each $m$, or the Cartesian reciprocal vectors, so round-tripping is unambiguous.
2. Define and document a canonical $G$-vector ordering for TREXIO.
   For example, sort by $|\mathbf{G}|$ and break ties lexicographically by $(h,k,l)$, and require all writers/readers to enforce the same convention.

Which approach do you think better?

Defining the plane-wave basis and the k-dependent number of PWs

Another challenge concerns how TREXIO defines the plane-wave basis. In §3.1 Basis set (basis group), the only PW-related field appears to be e_cut. It may be clearer to define the PW explicitly via the set of G-vectors (e.g., their Miller indices or Cartesian components), as discussed above. Also, PWs have no angular component; the functions $e^{,i\mathbf{G}_m\cdot\mathbf{r}}$ are better regarded as atomic orbitals rather than basis functions.

Relatedly, because the cutoff is applied to $|\mathbf{k}+\mathbf{G}|$, the number of PWs generally depends on $\mathbf{k}$, which we also confirmed from QE output files. TREXIO's AO group (4.1 Atomic orbitals) expects a single num for the number of functions, which makes $k$-dependent PW counts awkward unless we either (a) pad the MO coefficient arrays (4.2 Molecular orbitals: coefficient and coefficient_im) to a common length, or (b) allow per $\mathbf{k}$ arrays (but (b) breaks the backward compatibility?)

We don't have a firm proposal yet and would appreciate feedback from you two.

Thank you in advance.

[scemama]

Dear @kousuke-nakano @q-posev and @selim-tunbul,
I have no experience in periodic calculations, so it is likely that what we propose in trexio is incomplete.
Can we organize a zoom meeting next week to talk about this?

[kousuke-nakano]

Dear @scemama, Yes, thank you! Let's organize a Zoom meeting.

[selim-tunbul]

@scemama , after some discussion me(@selim-tunbul) and @kousuke-nakano have identified two possible approaches:
(A) Introducing a new section (PW)
New variables or sections: A new section “PW” containing XXX

Advantages:

• Does not require modification of existing structures and should remain backwards compatible.
• Provides a clearer separation between localized-orbital conventions and plane -waves-based conventions
• Relatively easy to implement

Disadvantages:

• In the future, we may want to use one-electron integrals( ao_1e_int group) and two-electron integrals(ao_2e_int group), but these groups depend on ao.num, which would not be defined by using the new pw group.

(B) Introducing the miller index array into the existing section (basis set)
New variables or sections: Miller index into the section Basis set. Nothing else changes.

In this approach, we would use mo.num = number of bands and ao.num =[to be defined].

Advantages:

• Existing structures would immediately allow the use of one-electron integrals( ao_1e_int group) and two-electron integrals(ao_2e_int group)
• No need to add extra group
  Disadvantages:
• This method is noy yet fully developed and and appears to require a more complicated implementation by developers.
• May reduce clarity by users, as plane-wave-related information would reside in the molecular orbitals group.

Additional note

• I was also thinking that since someone previously added or suggested the section on plane-wave implementation, if it might be possible for that person to review our proposal. Our current issue arises from a project to write a converter specifically for Quantum ESPRESSO. However, any recommendations we make should remain generally applicable, not limited to this specific case, and input from others would be valuable.
• However, for the plan(B), we need to consider the different lengths of Miller indices for each k. How can we do this?

Plan A:

Basis group

Plane wave
A plane-wave basis function is

$$\chi_{m}(\mathbf{r})=\exp(-i\mathbf{G}_{m}\cdot \mathbf{r})$$

where $\mathbf{G}_{m}$ is a reciprocal-lattice vector is

$$\mathbf{G}_{m}=h_{m}\mathbf{g}_{a}+k_{m}\mathbf{g}_{b}+l_{m}\mathbf{g}_{c}.$$

The triplet $(h_{m},k_{m},l_{m})$ gives the integer Miller indices of the $m$-th plane wave for a particular crystal momentum $\mathbf{k}$.

For band $\mathbf{n}$ at k-point $\mathbf{k}$, the Bloch state is expanded as

$$\psi_{n\mathbf{k}}(\mathbf{r})=\sum_{m}C_{n,m}^{\mathbf{k}}\chi_{m}(\mathbf{r})$$

The expansion coefficients $C_{n,m}^{k}$ are stored as two separate real arrays: coefficient and coefficient_im respectively. The ordering of both the Miller indices and the coefficients must match exactly:

$$\text{mills} = \begin{bmatrix} \boxed{ \begin{array}{ccc} h_{1} & k_{1} & l_{1}\\ \vdots & \vdots & \vdots \\ h_{M_1} & k_{M_1} & l_{M_1} \end{array}} \text{for $k_1$} \\ \boxed{ \begin{array}{ccc} h_{M_1+1} & k_{M_1+1} & l_{M_1+1}\\ \vdots & \vdots & \vdots \\ h_{M_1+M_2} & k_{M_1+M_2} & l_{M_1+M_2} \end{array} } \text{for $k_2$} \\ \vdots \end{bmatrix}$$

All coefficients are flattened into a single 1-D array, ordered first by band, then by $\mathbf{k}$-point, and finally by $\mathbf{G}$-vector:

$$\text{Coefficient}= [\overbrace{ C_{1,1}^{k_{1}},\dots,C^{k_{1}}_{1,M_{1}},C_{1,1}^{k_{2}},\dots }^{ \text{all k's for band 1} },\overbrace{ C_{2,1}^{k_{1}} }^{ \text{band 2} },\dots] $$

Variable	type	Row-major dimensions	Description	Implemented	TODO
type	str		Type of basis set: "PW"
total_num_pw	dim		Total number of plane waves		add
prim_num_kpoint	int	[pbc.k_point_num]	Number of plane waves per k-point		add
mills	int	[pbc.k_point_num $\times$ basis.prim_num, 3 ]	Miller indices of plane waves		add
e_cut	float		Energy cut-off for plane-wave		review

PW group

Variable	type	Row-major dimensions	Description	Implemented	TODO
coefficient	float	[pbc.k_point_num $\times$ basis.total_num_pw]	Real parts of $C_{n,m}^{k}$		add
coefficient_im	float	[pbc.k_point_num $\times$ basis.total_num_pw]	Imaginary parts of $C_{n,m}^{k}$		add

[kousuke-nakano]

Dear @selim-tunbul

First of all, thank you for the note! This is quite helpful for this issue! The PW basis set is:

Overview of the Plane wave basis

A plane-wave basis function is

$$\chi_{m}(\mathbf{r})=\exp(-i\mathbf{G}_{m}\cdot \mathbf{r})$$

where $\mathbf{G}_{m}$ is a reciprocal-lattice vector is

$$\mathbf{G}_{m}=h_{m}\mathbf{g}_{a}+k_{m}\mathbf{g}_{b}+l_{m}\mathbf{g}_{c}.$$

The triplet $(h_{m},k_{m},l_{m})$ gives the integer Miller indices of the $m$-th plane wave for a particular crystal momentum $\mathbf{k}$.

For band $\mathbf{n}$ at k-point $\mathbf{k}$, the Bloch state is expanded as

$$\psi_{n\mathbf{k}}(\mathbf{r})=\sum_{m}C_{n,m}^{\mathbf{k}}\chi_{m}(\mathbf{r})$$

The expansion coefficients $C_{n,m}^{k}$ are stored as two separate real arrays: coefficient and coefficient_im respectively. The ordering of both the Miller indices and the coefficients must match exactly:

$$\text{mills} = \begin{bmatrix} \boxed{ \begin{array}{ccc} h_{1} & k_{1} & l_{1}\\ \vdots & \vdots & \vdots \\ h_{M_1} & k_{M_1} & l_{M_1} \end{array}} \text{for $k_1$} \\ \boxed{ \begin{array}{ccc} h_{M_1+1} & k_{M_1+1} & l_{M_1+1}\\ \vdots & \vdots & \vdots \\ h_{M_1+M_2} & k_{M_1+M_2} & l_{M_1+M_2} \end{array} } \text{for $k_2$} \\ \vdots \end{bmatrix}$$

All coefficients are flattened into a single 1-D array, ordered first by band, then by $\mathbf{k}$-point, and finally by $\mathbf{G}$-vector:

$$\text{Coefficient}= [\overbrace{ C_{1,1}^{k_{1}},\dots,C^{k_{1}}_{1,M_{1}},C_{1,1}^{k_{2}},\dots }^{ \text{all k's for band 1} },\overbrace{ C_{2,1}^{k_{1}} }^{ \text{band 2} },\dots] $$

A possible solution

I think I have come up with a straightforward idea to solve the issue we discussed.

Basis group

I do not think we have to introduce something into the basis group.

Variable	type	Row-major dimensions	Description	Implemented	TODO
type	str		Type of basis set: "PW"

AO group

One should set ao.num = total number of plane waves.

MO group

Then, we can introduce a variable millar index in the MO group.

Variable	type	Row-major dimensions	Description	Implemented	TODO
k_point	index	[mo.num]	For periodic calculations, the $k$ point to which each MO belongs
millar index	int	[mo.num, 3 ]	Miller indices of plane waves for all $k$ points		add

[scemama]

Thanks for your effort!
I think @kousuke-nakano's proposition seems reasonable, with very minimal changes.