Skip to content

Conversation

mikaelkuisma1
Copy link

@mikaelkuisma1 mikaelkuisma1 commented Nov 6, 2022

I run the layer group detection to all 15000 2D-materials in C2DB. This resulted in 22 Bus errors or segfaults.

Here is a simple script modified from example.c to reproduce the problem.

gcc runlayers.c -I../include -L../_build -O0 -lsymspg -fsanitize=address -static-libasan -ggdb && LD_LIBRARY_PATH=pwd/../_build ./a.out

#include "spglib.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static void test_layer_spg_get_dataset(void);
static void show_layer_spg_dataset(double lattice[3][3],
                                   const double origin_shift[3],
                                   double position[][3],
                                   const int num_atom,
                                   const int types[],
                                   const int aperiodic_axis,
                                   const double symprec);

int main(void)
{
  test_layer_spg_get_dataset();
  return 0;
}

static void test_layer_spg_get_dataset(void)
{
  double symprec;
    double origin_shift[3] = { 0.0, 0.0, 0.0};
    int num_atom = 6;

    //double lattice[3][3] = { { 3.6298283346442517, 4.414688373086613e-05, 4.931753153254027e-16 }, { 2.0578765260243217e-05, 6.3189524931635095, 3.0514353036695804e-21 }, { 0.0, 0.0, 16.2670050174248 } };
    double lattice[3][3] = { { 3.6298283346442517, 2.0578765260243217e-05, 0.0 }, { 4.414688373086613e-05, 6.3189524931635095, 0.0 }, { 4.931753153254027e-16, 3.0514353036695804e-21, 16.2670050174248 } };
     int types[] = { 49, 49, 16, 16, 16, 16 };

    double position[][3] = {
{ 0.5000000000011727, 1.9659104439318793e-12, 0.499999999990874 },
{ 3.918861871595479e-13, 0.5000000000025439, 0.499999999990813 },
{ 0.49999671789715555, 0.3331107557679678, 0.4034517768684944 },
{ 0.9999871413290102, 0.16687354536345317, 0.5966360513378329 },
{ 0.5000032821094088, 0.6668892442284696, 0.5965482231161988 },
{ 1.2858669915597997e-05, 0.8331264546340158, 0.4033639486465312 }
 };

    show_layer_spg_dataset(lattice, origin_shift, position, num_atom, types, 2, 0.1);
}

static void show_layer_spg_dataset(double lattice[3][3],
                                   const double origin_shift[3],
                                   double position[][3],
                                   const int num_atom,
                                   const int types[],
                                   const int aperiodic_axis,
                                   const double symprec)
{
  SpglibDataset *dataset;
  char ptsymbol[6];
  int pt_trans_mat[3][3];

  int i, j, size;
  const char *wl = "abcdefghijklmnopqrstuvwxyz";

  for ( i = 0; i < num_atom; i++ ) {
    for ( j = 0; j < 3; j++ ) {
      position[i][j] += origin_shift[j];
    }
  }

  dataset = spg_get_layer_dataset(lattice,
                                  position,
                                  types,
                                  num_atom,
                                  aperiodic_axis,
                                  symprec);

  printf("International: %s (%d)\n", dataset->international_symbol, dataset->spacegroup_number );
  printf("Hall symbol:   %s\n", dataset->hall_symbol );
  spg_get_pointgroup(ptsymbol,
                     pt_trans_mat,
                     dataset->rotations,
                     dataset->n_operations);
  printf("Point group:   %s\n", ptsymbol);
  printf("Transformation matrix:\n");
  for ( i = 0; i < 3; i++ ) {
    printf("%f %f %f\n",
           dataset->transformation_matrix[i][0],
           dataset->transformation_matrix[i][1],
           dataset->transformation_matrix[i][2]);
  }
  printf("Wyckoff letters:\n");
  for ( i = 0; i < dataset->n_atoms; i++ ) {
    printf("%c ", wl[dataset->wyckoffs[i]]);
  }
  printf("\n");
  printf("Equivalent atoms:\n");
  for (i = 0; i < dataset->n_atoms; i++) {
    printf("%d ", dataset->equivalent_atoms[i]);
  }
  printf("\n");

  for (i = 0; i < dataset->n_operations; i++) {
    printf("--- %d ---\n", i + 1);
    for (j = 0; j < 3; j++) {
      printf("%2d %2d %2d\n",
             dataset->rotations[i][j][0],
             dataset->rotations[i][j][1],
             dataset->rotations[i][j][2]);
    }
    printf("%f %f %f\n",
           dataset->translations[i][0],
           dataset->translations[i][1],
           dataset->translations[i][2]);
  }

  spg_free_dataset(dataset);

}

I traced the error down to cel_layer_is_overlap function, which typically had periodic axes 0, 1, until they were 2 1070703038, and the code crashed. It is quite miraculous that the code worked before, probably due to some stack alignment of previous calls, that the periodic axes happened to be 0 and 1.

This merge request allows calculation for those 22 systems, while keeping all other results the same. However, I am still not sure are the results correct.

There are some more bugs remaining still. Edit: all of these might be resolved, see comment below.

  1. After changing the periodic axis to first axis from last axis, one material (from all 15000) went from LG 10 to LG 18, otherwise data was reproduced.
    periodic axes are hard coded in some places, however, I do not know if this is related
hall_symbol.c:    int periodic_axes[2] = {0, 1};
site_symmetry.c:    int periodic_axes[2] = {0, 1};
site_symmetry.c:    int periodic_axes[2] = {0, 1};
site_symmetry.c:    int periodic_axis[2] = {0, 1};
  1. There are 191 materials from 15k, where the layer group maps to other spacegroup than is given by get_spacegroup.
    (Mapped according to this table) https://arxiv.org/pdf/1306.1280.pdf
    Here LG is layergroup, and SG is spacegroup, and mappedsg is mapping of layer group to spacegroup by the table.
    https://pastebin.com/sUqCtu4j

Most of them are predicted to be layergroup 35, according to article, spacegroup 35, but predicted to be spacegroup 38.

@mikaelkuisma1
Copy link
Author

mikaelkuisma1 commented Nov 6, 2022

Update:
If I change the tolerance to 0.05 from 0.1, the all axes permutations pass on the one system which went from when changing the aperiodic axis previously. Thus, maybe there are multiple choices on this system, due to large tolerance.

Using an online tool, if I cut spacegroup 38, I get layergroup 35. Thus, if I change (lg 35 -> sg 35, to lg 35-> sg 38), along with the tolerance update: all 15733 materials pass the sanity check.

http://dcwww.camd.dtu.dk/~kuisma/spacegroup38.png

With these changes, the data is here:
http://dcwww.camd.dtu.dk/~kuisma/layergroups_update_spacegroup2.txt

@atztogo
Copy link
Collaborator

atztogo commented Nov 7, 2022

Thanks @mikaelkuisma1 for your bug fix along the detailed explanation (helpful!), and systematic test. How did you run over C2DB 15000 data systematically?

Currently I am busy. I would like to have time to look at it next week or later. But I have not yet been very familiar with the layer group, and so I need also some basic study on it, too.

@site-g, can you have a look at it? Probably we need to refactor the code, too. @lan496, do you have any comment? (also magnetic layer group for the future)

@site-g
Copy link

site-g commented Nov 7, 2022

  1. The obvious bug seems obviously a bug. I'm very sorry for that.

  2. I think the hard coded periodic axes won't cause the change. The periodic axes are hard coded because the aperiodic axis has
    been aligned to c axis in function get_axes when search for layer group.

search_hall_number

As you mentioned tolerance would alter the result, cause of the change might be due to some of the atoms staying at some subtle positions. Could you check whether the **spacegroup** would change if you manually rotate the structure to different axis in the structure file and search their spacegroup? As I use almost the same algorithm to find layer group.
  1. I'm afraid the mapping in 1306.1280 might not be correct. SG 35 Cmm2 and SG 38 Cm2m(C2mm) are different group, neither are their layer group counterpart LG 26 cmm2 and LG 35 cm2m(c2mm).
  • In SG/LG Cmm2, the two-fold axis is perpendicular to the C-center face, while in SG/LG Cm2m, it is parallel to the plane.
  • Besides, in LG cmm2, the two-fold axis will along the aperiodic axis, while it is perpendicular to the aperiodic axis in cm2m.
  • You can also verify this by comparing their symmetry operations:

@lan496
Copy link
Member

lan496 commented Nov 7, 2022

(Just a comment) For testing, we may need another repository to check layer-group implementations. For example, I've created and checked magnetic-space-group implementations for MAGNDATA in a private repo. For refactoring, wrapper functions to hide comparison of cell->aperiodic_axis may be helpful, like inline int is_aperiodic_axis(const *Cell cell, const int axis).

@codecov-commenter
Copy link

Codecov Report

Base: 85.12% // Head: 85.12% // No change to project coverage 👍

Coverage data is based on head (6df1f54) compared to base (dd4b0cf).
Patch has no changes to coverable lines.

Additional details and impacted files
@@           Coverage Diff            @@
##           develop     #199   +/-   ##
========================================
  Coverage    85.12%   85.12%           
========================================
  Files           16       16           
  Lines         1318     1318           
========================================
  Hits          1122     1122           
  Misses         196      196           

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report at Codecov.
📢 Do you have feedback about the report comment? Let us know in this issue.

@site-g
Copy link

site-g commented Nov 7, 2022

That sounds very helpful for avoiding typo!

In fact, the layer-group implementation is far from being completed. There requires:

  • Systematic test (Almost done by @mikaelkuisma1)
  • Something in datasets to indicate group type (now negative hall number means layer group)
  • C-APIs (now only spg_get_layer_dataset)
  • Python, etc. interfaces
  • Documentation for layer group

And also some codes need to be refactored, but I have being very busy since the end of last year and contribute nothing to the project. Very sorry.

@atztogo
Copy link
Collaborator

atztogo commented Nov 7, 2022

@site-g, no need to be sorry. Just providing some comments is helpful. So can I merge this PR? Then we will discuss on the layer group in a new issue.

@site-g
Copy link

site-g commented Nov 7, 2022

I think we can, this is a typo.

@atztogo atztogo merged commit 876fb9d into spglib:develop Nov 7, 2022
@atztogo atztogo mentioned this pull request Nov 7, 2022
6 tasks
@atztogo
Copy link
Collaborator

atztogo commented Nov 7, 2022

Merged. Thanks @mikaelkuisma1.
I made an issue #200 to continue the discussion. Thanks @site-g and @lan496 for your valuable comments.

@LecrisUT LecrisUT added the bug label Mar 1, 2023
@LecrisUT LecrisUT changed the title Draft: Fix a bug in layer groups. Fix a bug in layer groups. Mar 1, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants