materialsproject / pymatgen

Python Materials Genomics (pymatgen) is a robust materials analysis code that defines classes for structures and molecules with support for many electronic structure codes. It powers the Materials Project.
https://pymatgen.org
Other
1.52k stars 868 forks source link

Slab generator doesn't find all surfaces, possible bug in symmetry routines #1294

Closed MLars closed 6 years ago

MLars commented 6 years ago

Dear community,

first of all I would like to thank you for the great work you put into the project. I started working with pymatgen a few months ago. I stumbled across the "bug" by accident and did quite some tests for verification.

The intention to use pymatgen was the availability of a slab generator for crystals in order to study surfaces of hexagonal crystals, namely LiNiO2, which is a common material for batteries. I am using VASP, so I attached a POSCAR file and also minimal script which is able to reproduce the bug.

Description of the bug Surfaces of crystals are characterized by their Miller index in the conventional standard cell. Some crystal surfaces and thus, also the corresponding Miller indices, are related to each other by symmetry operations. The provided script uses pymatgen to determine all distinct Miller indices of crystal surfaces up to a maximum index of 2. This is the result (just grep for "Miller" in the output and remove the Miller indices occurring twice): Miller index: (0, 0, 1) Miller index: (1, 0, 0) Miller index: (1, 0, 1) Miller index: (1, 1, 0) Miller index: (1, 1, 1) Miller index: (1, 0, 2) Miller index: (2, 0, 1) Miller index: (2, 1, 0) Miller index: (2, 1, 1) Miller index: (2, 2, 1) Miller index: (2, 1, 2) Miller index: (2, -1, 2)

I was particularly interested in the plane with Miller index (1,0,2) and (0,1,2). By cleaving the crystal along these planes (use e.g. VESTA and the conventional cell!) it is clearly visible, that both planes are completely different. (1,0,2) is kind of a vicinal surface, whereas (0,1,2) is terminated with either a (Li-Ni) plane or an Oxygen plane. (1,0,2) surface in the top left (green plane): poscar_std_102

(0,1,2) surface in the top left (green plane): poscar_std_012

So my basic question is, why the (0,1,2) surface is not found by the SlabGenerator?

Possible root cause of the faulty behavior: I did some more investigations and found out, that there are some symmetry operations that do not map the crystal onto itself. With the commands (to be used in the script attached to this report) sga=SpacegroupAnalyzer(pin.structure) 'conv=sga.get_conventional_standard_structure()' symm_ops = SpacegroupAnalyzer(conv).get_symmetry_operations() it is possible to analyze the symmetry operations. In particular the symmetry operation in real space ((1,-1,0),(1,0,0),(0,0,1)) (with some finite translation tau), which is in reciprocal space ((0,-1,0),(1,1,0),(0,0,1)) would map the Miller index (1,0,2) to (0,1,2). However, as stated above, these planes are not equivalent. So it might be possible that I (and maybe also the programmers of the SlabGenerator) do not fully understand how symmetry works or is implemented, or there is some kind of bug in the symmetry routines of pymatgen or, in the worst case, even in spglib.

POSCAR file: POSCAR_LiNiO2.txt

script.py to generate slabs from crystals (please remove .txt): surface_minimal.py.txt

I would be very thankful for any feedback.

Best regards, Lars

shyuep commented 6 years ago

Thank you for the comprehensive report. However, I don't think you are doing the slabs right in VESTA. First of all, you can logically reason (without any code) that it is an impossibility that (102) is not equivalent to (012) in hexagonal systems. By definition, hexagonal systems have a 3/6-fold rotational symmetry. So anything that involves permutation of just the first two miller indices must result in symmetrically equivalent planes. In your case, you can simply think of it as the a direction in a hexagonal crystal is the same as the b direction.

I took the liberty of generating the (102) and (012) planes in the crystal, setting the d = 1. Download and open https://www.dropbox.com/s/wp1e0ii5zc8iskb/LNO.vesta?dl=0

You can quite clearly see that the two planes are equivalent.

MLars commented 6 years ago

Dear shyuep,

thank you for your time and your kind reply. First of all, I am aware of the fact that in hexagonal systems the 4-index notation should be valid and I would conclude that (102) and (012) are equivalent. The exact crystal system based on the spacegroup analysis of pymatgen is trigonal. But also for trigonal systems the 4-index notation should be valid (please correct me if I am wrong). This is also one of the reasons why I was surprised by the differences I've seen between the (102) and (012) surface.

I think I prepared the surfaces in VESTA correctly, I just used a slightly higher "d" and I repeated the cell several times in all directions to get a better view on the surfaces. I also took a look into your VESTA file (together with a colleague) and unfortunately we do not agree that the equivalence is " quite clear". To us both planes are not equally occupied with atoms, even if the plane is shifted along its normal direction I can't clearly see it.

I would really appreciate if you could go more into detail, how you concluded to "equivalent" planes based on the visual inspection in VESTA.

If you are experienced with pymatgen you could also generate some slabs in (102) and (012) direction and they look different.

In the meanwhile we also compared with another software that is used in our group and it also returns two different surfaces.

With best regards, Lars

shyuep commented 6 years ago

Sorry, you are right - the two planes are not equivalent. Let me ponder the problem a bit more.

shyuep commented 6 years ago

Lars, I did some basic investigations (note: I was the one who coded the slab generators). It seems to me that the hex setting of the rhom system does not preserve all symmetries. If you instead use the primitive rhom cell of LiNiO2 to generate the slabs, you have 35 slabs within a max miller index of 2 instead of 31 generated for the hexagonal non-primitive cell. If you map the (102) and (012) miller indices to the rhom miller indices, they are quite clearly non-equivalent in the rhom setting (which should follow similar permutation symmetries as the cubic system). I need to think about this more deeply before giving you a definite answer.

For now, can you try using the prim cell to generate the slabs and see if it reproduces the two non-equivalent slabs?

MLars commented 6 years ago

Dear Shyue Ping Ong,

is there a simple way in pymatgen to find the corresponding Miller indices in the primitive cell when the indices are known in the conventional cell? Furthermore, I am not able to reproduce the "35 slabs within a max miller index of 2", since it also depends on the settings of the slab generator, like "symmetrize" etc.
If I would know which Miller index is (012) but in the rhombohedral system it would be easier to check.

Best regards, Lars

shyuep commented 6 years ago

Yes. See solution to Q5 of my crystallography exam. You can perform a coordinate transformation to map Miller indices.

Spring 2017 Final Solutions.pdf

shyuep commented 6 years ago

Just a note that (012) in hex should map to (110) in rhom. (102) in hex should map to (411).

MLars commented 6 years ago

Today I confirmed that with the primitive cell as a starting point both surfaces are reproduced.

wenhaosun commented 6 years ago

Hi Lars,

You may also be interested in this recent publication.

https://www.sciencedirect.com/science/article/pii/S0039602817307537

There is a pymatgen implementation in the SI, but I haven't merged it into pymatgen yet.

Best, Wenhao Sun

On Thu, Oct 4, 2018 at 3:54 AM MLars notifications@github.com wrote:

Today I confirmed that with the primitive cell as a starting point both surfaces are reproduced.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/materialsproject/pymatgen/issues/1294#issuecomment-426974205, or mute the thread https://github.com/notifications/unsubscribe-auth/ABvGo3-Fr4y3F4VoUTF-zPjhEhTQmwhFks5uhejkgaJpZM4XD2FF .

-- Ceder Research Group http://ceder.berkeley.edu/ | Postdoctoral Fellow Lawrence Berkeley National Laboratory | Materials Sciences Division 1 Cyclotron Rd | Bldg 33-338F | Berkeley, CA 94704

MLars commented 6 years ago

Dear wenhaosun,

I am also aware of this publication and I also sent you an email regarding the stability of the program. I almost forgot about it, because I believe I fixed it for the time being, at least the results seemed reasonable.

Best regards, Lars

MLars commented 6 years ago

Dear Shyue Ping Ong,

is somebody already working on a fix for the bug or is there something I can/should do?

Best regards, Lars

shyuep commented 6 years ago

Yes we have some idea of a solution. @richardtran415 is working on it.

shyuep commented 6 years ago

This is fixed.

MLars commented 6 years ago

Thanks a lot. I will test it tomorrow.