fix #377 (symmetry): fix matrix pattern for layer groups

[site-g]

Under the convention of a-b axes defining the 2D lattice, this commit constrains the linear part $\boldsymbol{W}_\mathrm{layer}$ of layer group symmetry operations under the primitive cell basis after Delaunay reduction (and standardization) to

$$\begin{pmatrix} W_{11} & W_{12} & 0\\W_{21} & W_{22} & 0\\0 & 0 & \pm 1 \end{pmatrix} \tag{eq1}$$

Here, $W_{31} = 0$ because $\mathbf{a}'=W_{11}\ \mathbf{a}+W_{12}\ \mathbf{b} \pm \mathbf{c}$ will break the 2D lattice plane. By the same logic, $W_{32} = 0$.

To prove that $W_{13}$ and $W_{23} = 0$, we might need to enumerate all possibilities. It has been shown that, for 3D lattice, after Delaunay reduction, the shape of the primitive cell will belong to one of the 24 Symmetrische Sorten, and the transformation between them and the conventional cells are tabulated in Table 3.1.2.3 of ITA. For layer systems, the number of Symmetrische Sorten is smaller because the lattice plane is fixed. So, the primitive cell at this stage should be very close to the final standard primitive cell. And after selecting the three shortest vectors among the Delaunay set, the three vectors should be vertical to each other if it is still a primitive cell.

At this step, for layer groups with hexagonal, tetragonal and orthorhombic lattice systems, the basis vector $\mathbf{c}$ of the primitive cell should be vertical to the $\mathbf{ab}$ plane, so it is natural that $\boldsymbol{W}$s under these basis vectors have the form of (eq1) as $\mathbf{c}$ should keeps to be vertical after the symmetry operations . For the triclinic system, identity and inversion always satisfies the form (eq1).

Monoclinic is a little complicate. I think there are two allowed cases. First, $\mathbf{c}$ be vertical, and there is no problem. Second, $\mathbf{c}$ is inclined but being vertical to a 2-fold rotation whose rotation axis is in the $\mathbf{ab}$ plane. Such a rotation will change $\mathbf{c}$ to $-\mathbf{c}$ (Or, $\mathbf{c}$ is inclined but being vertical to the normal of a mirror plane whose normal is in the $\mathbf{ab}$ plane. Such a mirror reflection will leave $\mathbf{c}$ invariant), so (eq1) is still keeped.

What I showed in #377 is a cases that is not allowed. In that case, the in-plane component of $\mathbf{c}$ coincidentally equals $-\frac{1}{2}\mathbf{b}$, so is able to form a A-centred supercell and possess symmtry operations that are not allowed for layer groups. These illegal symmtry operations have non-zero $W_{13}, W_{23}$.

I think I have exhausted all possibilities, so the candidate operations should have the form (eq1).

[lan496]

@site-g Can you leave a short comment in the code to clarify we only care about rotation matrices whose non-digonal elements for the aperiodic axis are zero?
Then, can you run pre-commit?

pre-commit install
pre-commit run --all-files

[site-g]

I have added some comments. And the commit has passed all tests of pre-commit.

Besides, sorry for the force push in pull request #378. I was trying to create this pull request, but github cannot select commits and the new commit was automatically added to that pull request. Then I force push to drop the new commit.

[lan496]

Besides, sorry for the force push in pull request #378. I was trying to create this pull request, but github cannot select commits and the new commit was automatically added to that pull request. Then I force push to drop the new commit.

No problem! Git is difficult...

[site-g]

I have refactored the code. I used the very long if-clauses before because I thought we have known what are allowed and what are not allowed, so we can hardcode it rather than test it at runtime.

[codecov-commenter]

Codecov Report

All modified and coverable lines are covered by tests ✅

Comparison is base (8f7f01c) 83.80% compared to head (10f2bdc) 83.86%.
Report is 10 commits behind head on develop.

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #379      +/-   ##
===========================================
+ Coverage    83.80%   83.86%   +0.06%     
===========================================
  Files           24       24              
  Lines         8167     8180      +13     
===========================================
+ Hits          6844     6860      +16     
+ Misses        1323     1320       -3     

Flag	Coverage Δ
c_api	77.59% <100.00%> (+0.41%)	⬆️
fortran_api	56.15% <7.69%> (-0.05%)	⬇️
python_api	80.38% <38.46%> (-0.09%)	⬇️
unit_tests	1.24% <0.00%> (-0.01%)	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[LecrisUT]

About coverage, codecov is detecting that the other cases are not covered. Were they supposed to be in the python functional tests or ctest ones? Also specifically for these kind of logic, codecov only captures line coverage, while sonarcloud captures the cases inside the ifs as well¹. I will probably add this tool as well once I get to #364 #363.

For this PR though, if @lan496 knows where the relevant tests are, do you think you can cover at least the other aperiodic cases?

Footnotes

1. https://community.sonarsource.com/t/code-coverage-percentage-is-different-than-what-i-get-in-codecov/6527 ↩

[lan496]

The other cases are not covered because the unit test in https://github.com/spglib/spglib/blob/develop/test/functional/c/test_symmetry.cpp only uses aperiodic_axis=2, and there should be the place to test this PR.
I think just looping aperiodic_axis from 0 to 2 in test_get_lattice_symmetry_layer will cover the remaining lines.

@site-g Can you add the above tests?

[atztogo]

@lan496, I want to drop supporting aperiodic_axis = 0 or 1.

[lan496]

@atztogo Do you mean to drop aperiodic_axis from user inputs? Or even dropping aperiodic_axis from Cell struct?
If we can drop Cell->aperiodic_axis by imposing aperiodic_axis=2, I prefer it in a maintenance view.

[atztogo]

We haven't written the user documentation of layer group, I can say it is not released for general users.

What I wanted to mean was we don't need to spend time for writing tests for aperiodic_axis != 2 if all of you agree for dropping cases of aperiodic_axis != 2. We don't need to remove aperiodic_axis from codes right now, but we can just disable the cases of aperiodic_axis != 2 in any means. We can really remove it in the future versions at some big refactoring or in the process of the replacement of C code by some other language.

[site-g]

These codes were written when I was not familiar enough with the entire algorithm. Now I have realized that even if we still want to handle the aperiodic axis != 2 cases, we can transform it to be aperiodic axis == 2 at the very begin, and transform it back in the end. Then we only need to consider about aperiodic axis == 2 case in the middle the code.

So, I agree with Togo that we may not need to consider too much about the aperiodic axis != 2 cases. Anymore, I modified the test to cover the three different aperiodic axis choices.

[atztogo]

@site-g thanks for your work.

For monoclinic, we have the mp and op Bravais lattices. In mp, c is taken perpendicular to the lattice plane. In op, c is inclined generally. c is taken the shortest or second shortest for angle $\alpha$ obtuse. If the primary axis is 2 fold rotation axis, why Delaunay reduction is OK for finding c? Or at this step, you don't consider standardization yet?

[site-g]

@atztogo Because in the op case, Delaunay reduction makes sure that the bottom face being rectangular, i.e. $\mathbf{c}$ is inclined, but it is orthogonal to the 2 fold rotation axis (because it is the shortest non-lattice vector). So, after the rotation, $\mathbf{c}$ become $-\mathbf{c}$, and $W_{13}, W_{23}=0$.

[site-g]

However, I am sorry that (eq1) is still wrong. What we need to worry about is not op but oc and other conventional cell with A/B/body centring. In these case, $\mathbf{c}$ can be non-orthogonal with the 2 fold rotation, then $W_{13}, W_{23}$ can be nonzero.

So, I think we cannot solve #377 by simply changing $\boldsymbol{W}$. I made a new commit that treats the inclined $\mathbf{c}$ specially. It is based on the assertion that if $c$ is not vertical after Delaunay reduction, then the input structure can at most have a point group of 2/m. Then invalid operations are dropped.

[atztogo]

Thanks @site-g for your thought. I didn't think so deeply, and my question was simple. The job of symmetry.c is finding coset representatives. My guess was that (eq1) is OK. But you mention about standardization, which I thought it may not be achievable in symmetry.c. Therefore, my question was

1. How (or in which part of code) should the standardization be achieved? This is mainly the job out of symmetry.c.
2. Can (eq1) be enough for finding coset representatives? Or even if redundant, can we use, for example, the equation at the top of Bug: linear part W of layer group symmetry operations #377 safely?

In fact, my question is not well defined because we don't define how the job is divided into symmetry.c and out of symmetry.c. This is not only determined logically, but is also done by the design of the code.

[site-g]

1. Sorry, in fact I'm also not very sure what the "standardization" is. In Bug: linear part W of layer group symmetry operations #377, @atztogo you asked a question

3. Do you distinguish reduction and standardization?

Now I found that I may not understand that question correctly. Could you give more explanation about why you mentioned standardization there?

[atztogo]

Thank you for asking. Standardization that I mean is what is written here, https://spglib.readthedocs.io/en/latest/definition.html#spglib-conventions-of-standardized-unit-cell . This is not a proper terminology in crystallography, but loosely defined in spglib.

Do you distinguish reduction and standardization?

Yes. Delaunay reduction is sensitive to small numerical distortion, e.g., the result can depend on different hardware or perhaps compilers. The standardization to the conventional unit cell is more robust, at least for making the lattice being a Bravais lattice.

[site-g]

Then there is no standardization at all. We cannot standardize the cell before we know its lattice system, but that needs its space group.

What I mean is that Delaunay reduction has the effect to make the cell more close to a standard primitive cell. At least, if there is a aperiodic vector that is vertical to the lattice plane, then Delaunay reduction should choose it because vertical vector is always the shortest vector. I think if numerical distortion is big enough to affect these, then we should treat the aperiodic vector as inclined, i.e. the input can at most have a monoclinic symmetry.

[atztogo]

I see. This PR looks good to me.

[site-g]

2. Can (eq1) be enough for finding coset representatives? Or even if redundant, can we use, for example, the equation at the top of Bug: linear part W of layer group symmetry operations #377 safely?

I would like to answer this question and point out that this commit is still not satisfactory enough.

First, (eq1) will miss some symmtry operations. While

$$\boldsymbol{W} = \begin{pmatrix} W_{11} & W_{12} & W_{13} \\\ W_{21} & W_{22} & W_{23} \\\ 0 & 0 & \pm 1 \end{pmatrix} \tag{eq2}$$

will find illegal symmtry operations. So, none of them can be used safely. Here is an example (lattice vectors are column vectors)

  double lattice[3][3] = {{ 2.0, 0.0, -1.0},
                          { 0.0, 2.0,  0.0},
                          { 0.0, 0.0,  2.0}};
  double origin_shift[3] = {0.0, 0.0, 0.0};
  double position[][3] = {{ 0.0, 0.0, 0.0 },
                          { 0.5, 0.5, 0.0 }};
  int types[] = {1, 1};
  int num_atom = 2;

After Delaunay reduction, the primitive cell is

$$(\mathbf{a}_\mathrm{p},\mathbf{b}_\mathrm{p},\mathbf{c}_\mathrm{p})=\begin{pmatrix} -1.0 & -1.0 & 1.0\\-1.0 & 1.0 & 0.0\\0.0 & 0.0 & -2.0\end{pmatrix} \tag{eq3}$$

Fig. 1 A monolayer cut from an AB stacking square lattice material. The gray atoms are actual atoms. White atoms are virtual atoms because $\mathbf{c}$ does not have translation symmetry. $O,\ x$ and $y$ indicate the Cartesian coordinate system. The black cell is the Delaunay primitive cell. $\mathbf{a},\ \mathbf{b}$ and $\mathbf{c}$ indicate the basis vectors of the primitive cell. The blue cell is the input cell. Rotate the blue cell by 90° along $z$ get the red cell, but both of them generate same lattice. The blue and red plane are mirror plane that can be represented by integer matrix in the blue and red cell, repectively.

The 3D space group of the structure is I 4/m m m. So, when we search symmetry operations using (eq2) under the primitive basis vectors (eq3), we will find 4-fold rotation like

$$\begin{pmatrix} 0 & -1 & 1\\1 & 0 & 0\\0 & 0 & 1\end{pmatrix}\tag{eq4}$$

whose rotation axis is $(1, 1, 2)$, i.e. we have to double $\mathbf{c}_\mathrm{p}$ to get an axis vertical to the $\mathbf{a}\mathbf{b}$ plane, which is wrong.

If we search symmetry operations using (eq1), we will loss symmetry like

$$\begin{pmatrix} 0 & -1 & 1\\-1 & 0 & 1\\0 & 0 & 1\end{pmatrix}\tag{eq5}$$

which is a mirror reflection with normal $(1, 1, 0)$ along the $x$ direction (red plane in the figure).

So, neither of them can give correct results.

[site-g]

In the latest commit b173694 of this pull request. I fall back to use (eq2) because (eq1) misses symmetry operations. I first try to fix the problem of (eq2) by check the rotation axis direction of the found operations and drop the operation if its rotation axis cannot be represented by $(m, n, l)$, where $m, n \in \mathbb{N}$, $l = 0, \pm 1$. I also check that there are at least two linear independent axis that are vertical to the rotation axis and can also be represent by the format of $(m, n, l)$, otherwise I also drop the rotation. Using this procedure, I successfully dropped the illegal 4-fold operations.

However, the symmetry operations remained are: $\pm 1, 2_x, 2_y, m_x$ and $m_y$. Obviously, they does not form a group because $m_x \circ m_y =m_z$. $m_z$ has been dropped because it needs a vertical $\mathbf{c}$ axis. As you can see in Fig. 1, the two vertical mirror can only exist in two different cells (red and blue cell).

So, what I do now to fix this problem is to assert that the structure only has a monoclinic symmetry. Then I fix the unique axis to be $x$ (or $y$) and drop the symmetry along $y$ (or $x$). But I do not think this is a very good solution.

[site-g]

I think we may need to reconsider #323, the example I showed above shows that if the direction of $\mathbf{c}$ affects the symmetry we found. The found symmetry operations may not form a group. I think this is a serious problem, especially it happens at the early stage of the algorithm.

Now, I drop some operations to make it a group, but this is not very reasonable. What we found is a subset of the total symmetry, and it has the chance to miss some valid symmetry operations.

I think we may need to verticalize the $c$ axis at the start stage, and transform it back in the end. Then drop symmetry operations that cannot be represented by an integer matrix.

[atztogo]

I think we may need to verticalize the c axis at the start stage, and transform it back in the end. Then drop symmetry operations that cannot be represented by an integer matrix.

Please do not run this direction and please take time.

[site-g]

I see. That maybe another pull request in some big refactoring, and need to be done very carefully.

Then about this pull request. In summary, after this fix, we will search operations using (eq2), and drop invalid operations until the remaining operatons form a group. When the aperiodic axis of the Delaunay primitive cell is vertical to the 2D lattice plane, (eq2) will reduce to (eq1) and all the symmetry operations will be found. When the aperiodic axis is inclined, we will find a triclinic or monoclinic layer group. It may not contain all the symmetry operations. But at least it should not throw an error or fall into endless loops.

[atztogo]

@site-g, it is difficult for me to accept (eq2). I think the implementation strategy you explain above is not simple, and will be difficult to maintain even if it is correct.

Layer group is a stabilizer subgroup of space group. So the same set of the candidate W matrices as space group can be used as the first step, and the operations that preserve the lattice plane among possible W matrices as the space group will be left for the layer group. This algorithm is probably simple and sounds.

[site-g]

@atztogo I agree the stabilizer story is a simple and sounds theoretical picture. But it is not enough for practical use. It only tell us we should find operations that preserve the 2D lattice plane, but it did not tell us how to treat the non-lattice vector. And the commits I pushed here is handling the non-lattice vector, i.e. the aperiodic axis.

There are many theoretical pictures. The stabilizer picture is one of them. All the pictures should be equivalent and be statisfied by the practical code implementation, otherwise our code is totally wrong! I think the pictures are just explanation of the code, so we need to concern it when we write the manuscript or documentation. We should choose a picture that is easy to understand. But when we write the real code, I think the picture is not a problem. As long as the code is correct, all the pictures should be statisfied automatically. We need to concern more about details.

[site-g]

The commit looks complicate because I am trying to dealing with an inclined aperiodic axis. If you think it is difficult to maintain. Shall we follow the suggestion from @lan496 in #377?

I think we can constrain that the input aperiodic axis is always perpendicular to other periodic axes because I cannot think of any practical cases when the aperiodic axis should be not perpendicular.

Then we can use (eq1) safely. What we need to do additionally is to tell the user that the current version of Spglib can only handle perpendicular aperiodic axis.

[atztogo]

@site-g, I don't discuss about pictures. I want to make formulae on paper and implemented code as close as possible for the maintainability. If you can't describe what is implemented in a way that other people can follow, the approach is probably not good. Then, we have to think about a better approach. A better approach is that mathematically well supported. Otherwise, many exceptional cases can arise because what we can imagine by our brain about crystal symmetry is limited. So the first important is mathematical support, then implementation. I feel now it is the time for us to think about the mathematics by taking time.

But I don't still understand why (eq1) can not cover inclined c-axis.

It only tell us we should find operations that preserve the 2D lattice plane, but it did not tell us how to treat the non-lattice vector.

What we should get in symmetry.c is the set of coset representatives, for which the set of "operations that preserve the 2D lattice plane" is fine. Here I assume 2D Delaunay reduction of a and b axes, and c axis is taken to be shortest. How to treat the non-lattice vector should be implemented in the following steps, which I meant a part of "standardization".

[site-g]

But I don't still understand why (eq1) can not cover inclined c-axis.

For example, for structure with layer group c m 1 1

Fig. 2 This is a top view, the gray hollow circles are lattice of the upper layer (do not exist). Delaunay cell is $(\mathbf{a}, \mathbf{b}, \mathbf{c})$, the grey line is the mirror plane

$$\boldsymbol{W} = \begin{pmatrix}0 & -1 & -1\\-1 & 0 & -1\\0 & 0 & -1\end{pmatrix}$$

Only under the basis vector $\mathbf{c}'$ will $W_{13}$ and $W_{23} = 0$.

[atztogo]

Thanks @site-g, I understand.

Then the best strategy is what I wrote above, i.e.,

1. Find space group operations.
2. Find stabilizer subgroup for which lattice plane is given by user as a and b axes.

[atztogo]

I write it again because I modified my phrase. The best strategy is what I wrote above, i.e.,

1. Find space group operations.
2. Find stabilizer subgroup for which stabilized plane is given by user as a and b axes.

[site-g]

Fig. 3 A square lattice with inclined $\mathbf{c}$.

2. Find stabilizer subgroup for which stabilized plane is given by user as a and b axes.

You can find that the 4-fold rotation stabilized the lattice plane. So this strategy cannot solve the problem.

[atztogo]

Thanks again for your comments with figure, @site-g.
The discussion became involved and probably here is not the place for finding solution. I will contact you by any other means.

[site-g]

Thanks for the discussion with @atztogo. After the discussion, we decide to use (eq1). It works well on a structure with vertical c. But may not find some of the symmetry operations if c is still inclined after Delaunay reduction. I add some comments in the code to warning the developers about it.

[LecrisUT]

Could you fix up the test a bit while you are touching the file? Also please rebase onto latest develop branch