qzhu2017 / PyXtal

A code to generate atomic structure with symmetry
MIT License
234 stars 59 forks source link

space group #171

Closed pcpatfire closed 2 years ago

pcpatfire commented 2 years ago

I am trying to use your library to calculate the XRD pattern of crystal structure. I am entering the basic cell using ase, as for instance:

NaCl_ase = ase.Atoms(
    symbols=['Na', 'Na', 'Na', 'Na', 'Cl', 'Cl', 'Cl', 'Cl'],
    cell=[
        (5.62,    0,    0),
        (   0, 5.62,    0),
        (   0,    0, 5.62)],
    scaled_positions=[
        (0.000000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.500000),
        (0.500000, 0.000000, 0.500000),
        (0.500000, 0.500000, 0.000000),
        (0.500000, 0.500000, 0.500000),
        (0.500000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.000000),
        (0.000000, 0.000000, 0.500000)],
    pbc=True)

Then I calculate the xrd pattern using:

        struc = pyxtal()
        struc.from_seed(NaCl_ase)
        xrd = struc.get_XRD(thetas=(0,90))

And everything works just fine and the xrd pattern looks correct. The problem is when I want to use the basic crystal cell defined as:

    NaCl_ase = ase.Atoms(
    symbols=['Na', 'Cl'],
    cell=[
        (5.62,    0,    0),
        (   0, 5.62,    0),
        (   0,    0, 5.62) ],
    scaled_positions=[
        (0.000000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.500000)],
    pbc=True)

In this case I should be able to add the correct spacegroup number (225) to struc, otherwise the xrd pattern is not correct.

Is it correct to use struc.group = Group(225)?

I tried, but the xrd calculation does not seem to take into account the space group number I was setting. Is there a way to make it working? Thanks

qzhu2017 commented 2 years ago

@pcpatfire I think this is because you did not build the right structure in ASE.

import ase

NaCl_ase = ase.Atoms(
    symbols=['Na', 'Na', 'Na', 'Na', 'Cl', 'Cl', 'Cl', 'Cl'],
    cell=[
        (5.62,    0,    0),
        (   0, 5.62,    0),
        (   0,    0, 5.62)],
    scaled_positions=[
        (0.000000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.500000),
        (0.500000, 0.000000, 0.500000),
        (0.500000, 0.500000, 0.000000),
        (0.500000, 0.500000, 0.500000),
        (0.500000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.000000),
        (0.000000, 0.000000, 0.500000)],
    pbc=True)
print(NaCl_ase)

NaCl_ase = ase.Atoms(
symbols=['Na', 'Cl'],
cell=[
    (5.62,    0,    0),
    (   0, 5.62,    0),
    (   0,    0, 5.62) ],
scaled_positions=[
    (0.000000, 0.000000, 0.000000),
    (0.000000, 0.500000, 0.500000)],
pbc=True)
print(NaCl_ase)

They correspond to two different structures. You probably need to do something when you call ase.Atoms so that ase will recognize the space group symmetry.

Atoms(symbols='Na4Cl4', pbc=True, cell=[5.62, 5.62, 5.62])
Atoms(symbols='NaCl', pbc=True, cell=[5.62, 5.62, 5.62])
pcpatfire commented 2 years ago

Of course the two structures are different. The first one, with all the atoms and their positions, has a P1 space group. The second one needs a space group to be explicitely specified to apply the proper symmetry operations, otherwise there is no way, with such a basic cell, to identify automatically the space group by ase (or any other software). One question is: is it possible in pyxtal to set explicitly the space group by using the Groupclass? Another question is: is this space group used by pyxtal to expand the cell using the symmetry operations before calculating the xrd pattern?

qzhu2017 commented 2 years ago

Well, PyXtal is not really designed for that.

That said, below is a possible solution to work around:

>>> from pyxtal import pyxtal
>>> from pyxtal.lattice import Lattice
>>> c = pyxtal()
>>> l = Lattice.from_para(5.62, 5.62, 5.62, 90, 90, 90, l_type='cubic')
>>> c.from_random(3, 225, ['Na', 'Cl'], [4,4], lattice=l, sites=[["4a"], ["4b"]])
>>> c

------Crystal from random------
Dimension: 3
Composition: Cl4Na4
Group: Fm-3m (225)
triclinic lattice:   5.6200   5.6200   5.6200  90.0000  90.0000  90.0000
Wyckoff sites:
    Na @ [ 0.0000  0.0000  0.0000], WP [4a] Site [4/m-32/m]
    Cl @ [ 0.5000  0.5000  0.5000], WP [4b] Site [4/m-32/m]
>>> xrd=c.get_XRD(thetas=(1,90))

In future, I will probably change c.from_random to c.build

pcpatfire commented 2 years ago

Thanks for your answer. It's an interesting work around.

In my app I have a P1-like description of the crystal (atoms and their positions - I am using a diffy.structure Structureclass). By using this information plus the space group, is there a way to get the basic cell with its Wyckoff sites to apply your solution?

qzhu2017 commented 2 years ago

What do you mean by basic_cell? Do you mean the primitive cell?

If you input a P1 ASE structure, pyxtal will automatically detect the symmetry and standardize the cell according to the symmetry if you use .from_seed() function.

pcpatfire commented 2 years ago

Yes, I mean the primitive cell. I was already using this method too (see my first example). Is there a way to 're-set' the space group manually in pyxtal, after the .from_seed function call?

qzhu2017 commented 2 years ago

When you reset the space group from P1 to Fm-3m, you also increase the number of atoms by 4 times. This doesn't seem like a reasonable structure manipulation. Also, you maybe able to do this change for (0, 0, 0) from P1 to Fm-3m, since (0, 0, 0) is a well defined WP in Fm-3m. But it is not always the case for any space group (e.g, P212121). What are you going to handle such cases?

I don't understand why you need this function. Maybe you can explain a bit more about your motivation.

pcpatfire commented 2 years ago

If you start from the primitive cell when you reset the space group you don't multiply the atoms. My motivation is that in same cases the automatic symmetry detection may fail. The steps I have in mind are the following:

  1. get a P1 structure (diffpy)
  2. find the primitive cell (?)
  3. set the space group manually (?)
  4. calculate the xrd pattern (pyxtal)
qzhu2017 commented 2 years ago

Can you give a concrete example in which the symmetry detection does not work well? I think you can always play with the tolerance to guess a reasonable symmetry as the starting point. Assuming the initial guess is not far from the target solution, you may do the folllowing

  1. search for the supergroup symmetry if you think there are too many XRD peaks
  2. search for the subgroup symmetry structure if you want to see some splitting of PXRD peaks

Setting an arbitrary space group symmetry for a primitive cell can be very tricky! PyXtal provides the manipulation of structure according to sub/super-group symmetry. The code is probably not fully debugged yet. However, I am happy to help if you want.

pcpatfire commented 2 years ago

I have a real case in which the automatic detection seems not working.

I derived an ase P1 structure from the CIF file of a calcite (taken from COD database, COD #1010962). The CIF file contains the following definitions:

_symmetry_space_group_name_H-M   'R -3 c RS'
_cell_angle_alpha                46.1
_cell_angle_beta                 46.1
_cell_angle_gamma                46.1
_cell_length_a                   6.36
_cell_length_b                   6.36
_cell_length_c                   6.36

This is the corresponding ase Atomscall:

calcite_cod = Atoms(
                symbols=['Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 
                         'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 
                         'Ca', 'Ca', 'C', 'C', 'C', 'C', 'C', 'C', 'O', 
                         'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 
                         'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 
                         'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 
                         'O', 'O', 'O', 'O', 'O'],
                cell=[
                (4.18090, 1.87649, 4.41003),
                (0.0, 4.58270, 4.41003),
                (0.0, 0.0, 6.36),
                ],
                scaled_positions=[
                (0.25, 0.25, 0.25),
                (0.75, 0.0, 0.25),
                (0.0, 0.75, 0.25),
                (0.75, 0.75, 0.75),
                (0.25, 0.0, 0.75),
                (0.0, 0.25, 0.75),
                (0.9166666666666666, 0.5833333333333333, 0.5833333333333333),
                (0.41666666666666663, 0.3333333333333333, 0.5833333333333333),
                (0.6666666666666666, 0.08333333333333331, 0.5833333333333333),
                (0.41666666666666663, 0.08333333333333331, 0.08333333333333331),
                (0.9166666666666666, 0.3333333333333333, 0.08333333333333331),
                (0.6666666666666666, 0.5833333333333333, 0.08333333333333331),
                (0.5833333333333333, 0.9166666666666666, 0.9166666666666666),
                (0.08333333333333331, 0.6666666666666666, 0.9166666666666666),
                (0.3333333333333333, 0.41666666666666663, 0.9166666666666666),
                (0.08333333333333331, 0.41666666666666663, 0.41666666666666663),
                (0.5833333333333333, 0.6666666666666666, 0.41666666666666663),
                (0.3333333333333333, 0.9166666666666666, 0.41666666666666663),
                (0.0, 0.0, 0.0),
                (0.0, 0.0, 0.5),
                (0.6666666666666666, 0.3333333333333333, 0.3333333333333333),
                (0.6666666666666666, 0.3333333333333333, 0.8333333333333334),
                (0.3333333333333333, 0.6666666666666666, 0.6666666666666666),
                (0.3333333333333333, 0.6666666666666666, 0.16666666666666666),
                (0.25, 0.75, 0.0),
                (0.25, 0.5, 0.0),
                (0.5, 0.75, 0.0),
                (0.75, 0.25, 0.5),
                (0.5, 0.25, 0.5),
                (0.75, 0.5, 0.5),
                (0.75, 0.25, 0.0),
                (0.75, 0.5, 0.0),
                (0.5, 0.25, 0.0),
                (0.25, 0.75, 0.5),
                (0.5, 0.75, 0.5),
                (0.25, 0.5, 0.5),
                (0.9166666666666666, 0.08333333333333331, 0.3333333333333333),
                (0.9166666666666666, 0.8333333333333333, 0.3333333333333333),
                (0.16666666666666663, 0.08333333333333331, 0.3333333333333333),
                (0.41666666666666663, 0.5833333333333333, 0.8333333333333334),
                (0.16666666666666652, 0.5833333333333333, 0.8333333333333334),
                (0.41666666666666663, 0.8333333333333333, 0.8333333333333334),
                (0.41666666666666663, 0.5833333333333333, 0.3333333333333333),
                (0.41666666666666663, 0.8333333333333333, 0.3333333333333333),
                (0.16666666666666652, 0.5833333333333333, 0.3333333333333333),
                (0.9166666666666666, 0.08333333333333331, 0.8333333333333334),
                (0.16666666666666663, 0.08333333333333331, 0.8333333333333334),
                (0.9166666666666666, 0.8333333333333333, 0.8333333333333334),
                (0.5833333333333333, 0.41666666666666663, 0.6666666666666666),
                (0.5833333333333333, 0.16666666666666652, 0.6666666666666666),
                (0.8333333333333333, 0.41666666666666663, 0.6666666666666666),
                (0.08333333333333331, 0.9166666666666666, 0.16666666666666666),
                (0.8333333333333333, 0.9166666666666666, 0.16666666666666666),
                (0.08333333333333331, 0.16666666666666663, 0.16666666666666666),
                (0.08333333333333331, 0.9166666666666666, 0.6666666666666666),
                (0.08333333333333331, 0.16666666666666663, 0.6666666666666666),
                (0.8333333333333333, 0.9166666666666666, 0.6666666666666666),
                (0.5833333333333333, 0.41666666666666663, 0.16666666666666666),
                (0.8333333333333333, 0.41666666666666663, 0.16666666666666666),
                (0.5833333333333333, 0.16666666666666652, 0.16666666666666666),
                ],
                pbc=True)

When I call the .from_seed function it raises the following error:

Traceback (most recent call last):
  File "d:\work\CryStIE-devel\test_pyxtaldiff_calcitecod.py", line 212, in pmgdiff
    struc.from_seed(NaCl_ase)
  File "C:\Users\PC\AppData\Local\Programs\Python\Python39\lib\site-packages\pyxtal\__init__.py", line 337, in from_seed
    self._from_pymatgen(pmg_struc, tol, a_tol)
  File "C:\Users\PC\AppData\Local\Programs\Python\Python39\lib\site-packages\pyxtal\__init__.py", line 396, in _from_pymatgen
    raise RuntimeError("Cannot extract the right mapping from spglib")
RuntimeError: Cannot extract the right mapping from spglib

The problem is that the CIF does not contain the usual definition of the structure, as it can be found, for instance, in a calcite from the American Mineralogist database (code 0008879):

_cell_length_a 5.0492
_cell_length_b 5.0492
_cell_length_c 17.3430
_cell_angle_alpha 90
_cell_angle_beta 90
_cell_angle_gamma 120
_symmetry_space_group_name_H-M 'R -3 c'

Here the lattice is defined in the usual way, with alpha=beta=90 and gamma=120, and a=b. The corresponding P1 ase structure is:

calcite_ams = Atoms(symbols=['Ca', 'Ca', 'Ca', 'Ca', 'Ca', 'Ca', 
                'C', 'C', 'C', 'C', 'C', 'C', 'O', 
                'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 
                'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O', 'O'],
                cell=[
                (4.37273, -2.52459, 3.09174e-16),
                (0.0, 5.0492, 3.09174e-16),
                (0.0, 0.0, 17.343),
                ],
                scaled_positions=[
                (0.0, 0.0, 0.0),
                (0.6666666666666666, 0.3333333333333333, 0.3333333333333333),
                (0.3333333333333333, 0.6666666666666666, 0.6666666666666666),
                (0.0, 0.0, 0.5),
                (0.6666666666666666, 0.3333333333333333, 0.8333333333333334),
                (0.3333333333333333, 0.6666666666666666, 0.16666666666666666),
                (0.0, 0.0, 0.25),
                (0.6666666666666666, 0.3333333333333333, 0.5833333333333333),
                (0.3333333333333333, 0.6666666666666666, 0.9166666666666666),
                (0.0, 0.0, 0.75),
                (0.6666666666666666, 0.3333333333333333, 0.08333333333333348),
                (0.3333333333333333, 0.6666666666666666, 0.41666666666666663),
                (0.25642, 0.0, 0.25),
                (0.9230866666666666, 0.3333333333333333, 0.5833333333333333),
                (0.5897533333333334, 0.6666666666666666, 0.9166666666666666),
                (0.25642, 0.25642, 0.75),
                (0.9230866666666666, 0.5897533333333334, 0.08333333333333348),
                (0.5897533333333334, 0.9230866666666666, 0.41666666666666663),
                (0.0, 0.25642, 0.25),
                (0.6666666666666666, 0.5897533333333334, 0.5833333333333334),
                (0.3333333333333333, 0.9230866666666666, 0.9166666666666666),
                (0.74358, 0.0, 0.75),
                (0.41024666666666665, 0.3333333333333333, 0.08333333333333348),
                (0.07691333333333333, 0.6666666666666666, 0.41666666666666663),
                (0.74358, 0.74358, 0.25),
                (0.41024666666666665, 0.07691333333333333, 0.5833333333333334),
                (0.07691333333333333, 0.41024666666666665, 0.9166666666666666),
                (0.0, 0.74358, 0.75),
                (0.6666666666666666, 0.07691333333333333, 0.08333333333333348),
                (0.3333333333333333, 0.41024666666666665, 0.41666666666666663),
                ],
                pbc=True)

In this case the .from_seed function works fine. They have the same space group (#167), but the lattice is defined in a different way - this is possible for the calcite.

One note concerning your 'work around' to solve my problem. As you said you would change the .from_random function to .build. It would be great if it would be possible to create a structure using this kind of code:

>>> c = pyxtal()
>>> l = Lattice.from_para(5.62, 5.62, 5.62, 90, 90, 90, l_type='cubic')
>>> c.from_random(3, 225, ['Na', 'Cl'], [4,4], lattice=l, sites=[(0, 0, 0), (0.5, 0.5, 0.5)])

allowing to avoid the use of the Wyckoff positions, and making it possible to set no random numbers for the Wyckoff positions defined as (x,0,0), for instance, where the x value is not constant and it is assigned a random number by the .from_random function.

Last note is about the fact that in some cases, if the structure can have different origins and the space group is non-standard, it is changed to the standard one without any warning. I don't know if it possible to warn the user in some ways, but it would be useful.

Sorry for this long comment and thanks for your help.

qzhu2017 commented 2 years ago

@pcpatfire check #172 first. I will probably deal with your other issues later.

pcpatfire commented 2 years ago

Ok. Thanks

qzhu2017 commented 2 years ago

Why don't you just download the cif file and load it from pyxtal?

>>> from pyxtal import pyxtal
>>> c = pyxtal()
>>> c.from_seed('1010962.cif')
>>> c

------Crystal from Seed------
Dimension: 3
Composition: C6O18Ca6
Group: R-3c (167)
hexagonal lattice:   4.9803   4.9803  17.0187  90.0000  90.0000 120.0000
Wyckoff sites:
    Ca @ [ 0.0000  0.0000  0.0000], WP [6b] Site [-3..]
     C @ [ 0.0000  0.0000  0.2500], WP [6a] Site [32.]
     O @ [ 0.7500  0.0000  0.2500], WP [18e] Site [.2.]
>>> ase_obj = c.to_ase()
>>> ase_obj
Atoms(symbols='C6O18Ca6', pbc=True, cell=[[4.3130801187158605, -2.490157967577024, 3.049563984364482e-16], [0.0, 4.980315935154049, 3.049563984364482e-16], [0.0, 0.0, 17.01868853813807]])
pcpatfire commented 2 years ago

Yes. Of course you can load the downloaded cif. Buts as you can see the lattice is different, and if I plot the structure the cell does not correspond to the one defined in the original cif. 1010962-calcite-struc

I am creating an application that loads crystals in several different ways (not just from cif) and handles them as P1, because in general I could know anything else than the atoms' positions and lattice. Then I'm using the ase structure to calculate the xrd pattern using your.get_XRD function. But with some structures, as the calcite of the previous example, I get the error raised by .from_seed. I think the problem in that case is that a R-3c is associated to a romboidal cell with a=b=c and alpha=beta=gamma, instead of using the usual cell with a=b and two angles of 90° and one of 120°. I was using the calcite as a kind of blind test to see if I was able to handle any kind of correct description of a structure. Probably the problem is due to the spglib, as written in the error message. Unfortunately, from the point of view of my application it's a problem because, in principle, I should be able to calculate the xrd pattern for any structure, even if it is described in a "strange" way.

qzhu2017 commented 2 years ago

Based on your description, I think you can directly call our XRD module with the ASE-P1 structure. There is no need to transform it to the PyXtal object. PyXtal object is really designed for symmetry analysis.

from pyxtal.XRD import XRD
import ase

NaCl_ase = ase.Atoms(
    symbols=['Na', 'Na', 'Na', 'Na', 'Cl', 'Cl', 'Cl', 'Cl'],
    cell=[
        (5.62,    0,    0),
        (   0, 5.62,    0),
        (   0,    0, 5.62)],
    scaled_positions=[
        (0.000000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.500000),
        (0.500000, 0.000000, 0.500000),
        (0.500000, 0.500000, 0.000000),
        (0.500000, 0.500000, 0.500000),
        (0.500000, 0.000000, 0.000000),
        (0.000000, 0.500000, 0.000000),
        (0.000000, 0.000000, 0.500000)],
    pbc=True)

myxrd = XRD(NaCl_ase, thetas=[0,90])
print(myxrd.pxrd)
pcpatfire commented 2 years ago

Yes, it is a possibility, but for the calcite from COD I can't recreate the original group and cell.

If pyxtal is designed for symmetry analysis, why in the case of the calcite from the American Mineralogist database it can find the space group, while it raises an error with the one from COD? Maybe for the "strange" structures it is needed to take into account more information.

It is the same problem with the structures that have more than one possible description - one of my other issues - as for instance some structures with space group 62. In this case, if I load the olivine cif with COD number 1011090 (space group # 62, Pbnm):

_symmetry_space_group_name_Hall  '-P 2c 2ab'
_symmetry_space_group_name_H-M   'P b n m'
_cell_angle_alpha                90
_cell_angle_beta                 90
_cell_angle_gamma                90
_cell_formula_units_Z            4
_cell_length_a                   4.787
_cell_length_b                   10.086
_cell_length_c                   5.939
_cell_volume                     286.7
_cod_database_code               1011090

pyxtal transforms the space group into # 62, Pnma, which is the same group with modified lattice axis: a->c, c->b, b->a:

data_from_pyxtal

_symmetry_space_group_name_H-M 'Pnma'
_symmetry_Int_Tables_number                   62
_symmetry_cell_setting              orthorhombic
_cell_length_a           10.086000
_cell_length_b            5.939000
_cell_length_c            4.787000
_cell_angle_alpha        90.000000
_cell_angle_beta         90.000000
_cell_angle_gamma        90.000000
_cell_volume            286.744909

It is another possible description of the same crystal (the volume and the xrd pattern are the same). Why not ask the user the "preferred" one?

qzhu2017 commented 2 years ago

A crystal structure can have multiple representations. Rhombohedral structure can be expressed by either rhombohedral or hexagonal cell. In pyxtal, I choose the hexagonal setting, because I hate dealing with a cell that has a strange inclination angle.

For the 2nd structure, there exist 530 settings for the entire 230 space groups. http://cci.lbl.gov/sginfo/hall_symbols.html For each one of the 230 groups, Pyxtal only choose the standard setting. For space group 62, the standard setting is Pnma. If you swap the a,b,c axis, you will end up with different space group symbol. This is just an alternative representation. For most of applications, I think it is easier to deal with the standard setting.

It is possible to describe the structure with the customized setting, however, this will require some extra work.

I still don't understand why you care about the change of setting (if the structure does not change). If you really want this feature, I encourage you to implement it by yourself. PyXtal is an open source code. I am trying to make it easier to let newcomers join the development.

pcpatfire commented 2 years ago

Ok. Thanks