[Advanced] Structural optimization with symmetry¶
In the previous section, we learned the most primitive application of structural optimization. In this section, we will learn how to perform structural optimization while preserving the symmetry of the system. ASE has a method to handle space groups by using the constraints class FixSymmetry.
This section includes some advanced contents; those who wish to hurry to the next section may skip it and read it later.
BaTiO\(_3\) optimization¶
Let us look at a concrete example before going into the details of space groups and what we mean by “symmetry” in this section.
Here we will look at a structural optimization example of a material with a perovskite structure BaTiO\(_3\) (barium titanate).
BaTiO\(_3\) is an artificial mineral used as a dielectric material for ceramic laminated capacitors because of its high relative permittivity. It takes a tetragonal crystal system at room temperature. For more information, please refer to the following page.
The structure file is obtained from Materials Project.
The following cells of BaTiO\(_3\) with the tetragonal structure are \(a = b \neq c\), \(\alpha = \beta = \gamma =\) 90°. (\(a = b = c\), \(\alpha = \beta = \gamma =\) 90° case is called cubic)
[1]:
from ase.io import read
from pfcc_extras.visualize.view import view_ngl
bto_tetra = read("../input/BaTiO3_mp-5986_computed.cif")
print(f"BaTiO3 (Tetragonal): cell = {bto_tetra.cell.cellpar()}")
view_ngl(bto_tetra, representations=[])
BaTiO3 (Tetragonal): cell = [ 4.004457 4.004457 4.200636 90. 90. 90. ]
[2]:
view_ngl(bto_tetra * (4, 4, 4), representations=[])
The get_scaled_positions method can be used to obtain the relative coordinates of each atom in the unit cell, rather than the absolute coordinates. (That is, a coordinate scale such that [0, 0, 0] is the origin and [1, 1, 1] is the point at the opposite end of the unit cell.)
If you check the scaled_positions in this structure, you will see that we have taken a rounded number, such as 0 or 0.5, for the xy coordinates.
[3]:
bto_tetra.get_scaled_positions()
[3]:
array([[0.5 , 0.5 , 0.004136],
[0. , 0. , 0.524313],
[0. , 0.5 , 0.479349],
[0.5 , 0. , 0.479349],
[0. , 0. , 0.958854]])
Let’s try to optimize the structure in the following four ways for this structure,
With or without
FixSymmetryas described in this sectionExpCellFilter’shydrostatic_strainargument isTrueorFalse, which was described in the previous section.
[4]:
import pfp_api_client
from pfp_api_client.pfp.calculators.ase_calculator import ASECalculator
from pfp_api_client.pfp.estimator import Estimator, EstimatorCalcMode
estimator = Estimator(calc_mode=EstimatorCalcMode.CRYSTAL, model_version="v3.0.0")
calculator = ASECalculator(estimator)
Potential energy before structural optimization is
[5]:
bto_tetra.calc = calculator
Epot = bto_tetra.get_potential_energy()
print(f"Epot before opt: {Epot} eV")
Epot before opt: -31.54281016948673 eV
[6]:
from ase import Atoms
from ase.calculators.calculator import Calculator
from ase.constraints import ExpCellFilter
from ase.optimize import FIRE
from ase.constraints import FixSymmetry
def opt_with_symmetry(
atoms_in: Atoms,
calculator: Calculator,
fix_symmetry: bool = False,
hydrostatic_strain: bool = False,
) -> Atoms:
atoms = atoms_in.copy()
atoms.calc = calculator
if fix_symmetry:
atoms.set_constraint([FixSymmetry(atoms)])
ecf = ExpCellFilter(atoms, hydrostatic_strain=hydrostatic_strain)
opt = FIRE(ecf, logfile=None)
opt.run(fmax=0.005)
cell_diff = (atoms.cell.cellpar() / atoms_in.cell.cellpar() - 1.0) * 100
print("Optimized Cell :", atoms.cell.cellpar())
print("Optimized Cell diff (%):", cell_diff)
print("Scaled positions :\n", atoms.get_scaled_positions())
print(f"Epot after opt: {atoms.get_potential_energy()} eV")
return atoms
1. No FixSymmetry, hydrostatic_strain=False¶
The crystal axes a, b, c, \(\alpha, \beta, \gamma\) are all optimized independently, so you can see that the tetragonal symmetry is broken as the angles shift from 90°, although only slightly. The atomic coordinates are also moving freely.
[7]:
atoms1 = opt_with_symmetry(bto_tetra, calculator, fix_symmetry=False, hydrostatic_strain=False)
/tmp/ipykernel_1750/819195949.py:18: FutureWarning: Import ExpCellFilter from ase.filters
ecf = ExpCellFilter(atoms, hydrostatic_strain=hydrostatic_strain)
Optimized Cell : [ 3.999351 3.99935081 4.22621354 89.99999973 89.99999776 89.99999744]
Optimized Cell diff (%): [-1.27507905e-01 -1.27512654e-01 6.08896886e-01 -3.02345160e-07
-2.49438050e-06 -2.84332935e-06]
Scaled positions :
[[5.00000016e-01 5.00000010e-01 8.53036289e-03]
[1.59435257e-07 7.06001742e-08 5.26177155e-01]
[9.99999949e-01 4.99999954e-01 4.77548623e-01]
[4.99999884e-01 9.99999975e-01 4.77548659e-01]
[9.99999992e-01 9.99999991e-01 9.56196200e-01]]
Epot after opt: -31.543516660199565 eV
2. Apply FixSymmetry, hydrostatic_strain=False¶
The tetragonal symmetry of \(a = b \neq c\), \(\alpha = \beta = \gamma =\) 90° is well preserved. Also, we can see that the scaled_positions of the xy-coordinates also keep a number such as 0 or 0.5 (except for very small numerical errors). This means that the relationship between sites imposed by the symmetry P 4 m m described below is preserved. This is the desired result for this structural optimization.
[8]:
atoms2 = opt_with_symmetry(bto_tetra, calculator, fix_symmetry=True, hydrostatic_strain=False)
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['std_lattice']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['transformation_matrix']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['std_rotation_matrix']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['std_positions']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['mapping_to_primitive']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['std_mapping_to_primitive']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['rotations']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/home/jovyan/.py39/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['translations']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
warnings.warn(
/tmp/ipykernel_1750/819195949.py:18: FutureWarning: Import ExpCellFilter from ase.filters
ecf = ExpCellFilter(atoms, hydrostatic_strain=hydrostatic_strain)
Optimized Cell : [ 3.99935145 3.99935145 4.22621104 90. 90. 90. ]
Optimized Cell diff (%): [-0.12749674 -0.12749674 0.60883732 0. 0. 0. ]
Scaled positions :
[[5.00000000e-01 5.00000000e-01 8.53034739e-03]
[7.17464814e-43 3.00975353e-35 5.26176962e-01]
[0.00000000e+00 5.00000000e-01 4.77548635e-01]
[5.00000000e-01 0.00000000e+00 4.77548635e-01]
[0.00000000e+00 0.00000000e+00 9.56196419e-01]]
Epot after opt: -31.54352290834058 eV
3. No FixSymmetry, hydrostatic_strain=True¶
If hydrostatic_strain=True, the optimization preserves the ratio of lattice lengths and angles because it applies isotropic pressure. As you can see in the “Optimized Cell diff (%)” section, the angle is still 90°, but the a, b, and c axes are all optimized with the same ratio. As a result, the optimization of the structure extending in the direction of the c-axis has not progressed to a minimum value. (Compared to the cases 1 and 2, the structural optimization ends with a large value of
Epot.)
Since the atomic coordinates are free to move, the scaled_positions of the xy-coordinates have moved from a number such as 0 or 0.5.
[9]:
atoms3 = opt_with_symmetry(bto_tetra, calculator, fix_symmetry=False, hydrostatic_strain=True)
/tmp/ipykernel_1750/819195949.py:18: FutureWarning: Import ExpCellFilter from ase.filters
ecf = ExpCellFilter(atoms, hydrostatic_strain=hydrostatic_strain)
Optimized Cell : [ 4.00530335 4.00530335 4.20152382 90. 90. 90. ]
Optimized Cell diff (%): [0.0211353 0.0211353 0.0211353 0. 0. 0. ]
Scaled positions :
[[5.00000029e-01 5.00000017e-01 7.01999267e-03]
[1.64765465e-07 1.84306503e-07 5.24653378e-01]
[9.99999961e-01 4.99999863e-01 4.78119753e-01]
[4.99999878e-01 9.99999952e-01 4.78119744e-01]
[9.99999968e-01 9.99999985e-01 9.58088133e-01]]
Epot after opt: -31.54333787703485 eV
4. Apply FixSymmetry, hydrostatic_strain=True¶
For cells, the results are similar to those in 3. The constraint of maintaining the lattice length ratio and lattice angle as hydrostatic_strain=True is more stringent than those imposed by only FixSymmetry. For atomic coordinates, the FixSymmetry constraints keep the scaled_positions of the xy coordinates to 0 or 0.5.
[10]:
atoms4 = opt_with_symmetry(bto_tetra, calculator, fix_symmetry=True, hydrostatic_strain=True)
/tmp/ipykernel_1750/819195949.py:18: FutureWarning: Import ExpCellFilter from ase.filters
ecf = ExpCellFilter(atoms, hydrostatic_strain=hydrostatic_strain)
Optimized Cell : [ 4.00530331 4.00530331 4.20152377 90. 90. 90. ]
Optimized Cell diff (%): [0.0211343 0.0211343 0.0211343 0. 0. 0. ]
Scaled positions :
[[5.00000000e-01 5.00000000e-01 7.01994807e-03]
[1.76563607e-42 0.00000000e+00 5.24653369e-01]
[2.79559044e-43 5.00000000e-01 4.78119771e-01]
[5.00000000e-01 0.00000000e+00 4.78119771e-01]
[3.36311631e-43 0.00000000e+00 9.58088142e-01]]
Epot after opt: -31.543327059394755 eV
Summary¶
The results can be summarized as the below table.
Case |
|
|
Result |
|---|---|---|---|
1 |
Cell angles, etc., change slightly because symmetry is not taken into account. |
||
2 |
✓ |
The desired result. Structural optimization can be performed while maintaining the symmetry of a tetragonal crystal. The relationship between the atomic coordinate sites is also properly constrained. |
|
3 |
✓ |
Since all a, b, and c axes are optimized in the same ratio, the structural optimization in the c-axis direction is not sufficient. Atomic coordinates also move freely. |
|
4 |
✓ |
✓ |
Cell is similar to 3. The relationship of sites in atomic coordinates is subject to constraints. |
The symmetry is a “lattice vector & equivalent site relationship” (a=b, α=β=γ=90°, etc.), and using FixSymmetry allows structural optimization to be performed while preserving this symmetry. More precisely, structural optimization is performed only with displacements that preserve the space group to which the initial structure belongs. In the following, we will explain the definition of symmetry, i.e., the space group.
Space group (spacegroup)¶
It is known that there are a total of 230 types of patterns that can be taken by crystal structures that take regular pattern in three-dimensional space with periodic boundary conditions, which are called space groups. All crystals belong to one of the 230 space groups.
For a list of 230 species of space groups, you can check the following sites
The space group of atoms can be specified by using get_spacegroup in ASE.
[Note]
The library spglib is used internally in the get_spacegroup method. spglib is a library that determines the space group from unit cell and atomic coordinate information.
“
Spglib: a software library for crystal symmetry search”, Atsushi Togo and Isao Tanaka, https://arxiv.org/abs/1808.01590
Let’s determine the space group for this BaTiO3 tetragonal structure.
[11]:
from ase.spacegroup import get_spacegroup
sg = get_spacegroup(bto_tetra)
sg
[11]:
Spacegroup(99, setting=1)
The following information can be obtained from the Spacegroup class,
no: Number of space groupsymbol: Hermann-Mauguin notation. (Note that Schoenflies is the another famous notation.)lattice: Lattice type.Pis simple lattice,Ibody centered,Fface centered, etc.scaled_primitive_cell: Standardized primitive cell.reciprocal_cell: Reciprocal cell.
[12]:
type(sg)
[12]:
ase.spacegroup.spacegroup.Spacegroup
[13]:
print("no : ", sg.no)
print("symbol : ", sg.symbol)
print("lattice : ", sg.lattice)
print("scaled_primitive_cell: \n", sg.scaled_primitive_cell)
print("reciprocal_cell : \n", sg.reciprocal_cell)
no : 99
symbol : P 4 m m
lattice : P
scaled_primitive_cell:
[[1. 0. 0.]
[0. 1. 0.]
[0. 0. 1.]]
reciprocal_cell :
[[1 0 0]
[0 1 0]
[0 0 1]]
This BaTiO\(_3\) structure was determined to be the 99th P 4 m m space group.
The space group specifies the possible symmetry operations (such as translation, rotation, mirror image flips, and glide) that do not change its crystal structure. The details of the theory of space group, including how to read symbols, are out of the scope of this tutorial. Please refer to the references at the end of this document for those who wish to learn more.
In this tutorial, we will try to give an intuitive explanation of what space groups are by visualizing them. P 4 m m is a space group with 4-fold rotational symmetry (invariant to transformations that rotate it 90 degrees along the z-axis) and mirror-reversal symmetry with respect to the xz-plane and yz-plane.
For a given atomic configuration, the “equivalent sites” imposed by this symmetry can be calculated with the method equivalent_lattice_points. Let us calculate and display the equivalent sites for a point ([0.3, 0.1, 0.6]) as follows.
[14]:
sites, kinds = sg.equivalent_sites([[0.3, 0.1, 0.6]])
print("sites", sites)
p4mm_atoms = Atoms(symbols="C8", scaled_positions=sites, cell=[5, 5, 5], pbc=True)
view_ngl(p4mm_atoms)
sites [[0.3 0.1 0.6]
[0.7 0.9 0.6]
[0.9 0.3 0.6]
[0.1 0.7 0.6]
[0.3 0.9 0.6]
[0.7 0.1 0.6]
[0.9 0.7 0.6]
[0.1 0.3 0.6]]
Thus, we can see that a total of 8 points are listed as equivalent sites, including points obtained by rotating 90 degrees along the z-axis from [0.3, 0.1, 0.6] and by performing a mirror image flip from those points to the xz-plane/yz-plane.
Let us look again at the structure of BaTiO\(_3\).
[15]:
view_ngl(bto_tetra)
[16]:
from pfcc_extras.visualize.ase import view_ase_atoms
view_ase_atoms(bto_tetra, figsize=(4, 4), title="BaTiO3 tetragonal atoms", scale=100, rotation="0x,0y,0z")
As for the BaTiO\(_3\) structure in this case:
For Ba at index=0, Ti at index=1, and O at index=4, the only equivalent site is themselves since they move to the same location after rotation or mirror image inversion, except for a unit cell shift. The O structures at index 2 and 3 are found to be equivalent sites to each other due to 4-fold rotational symmetry.
[17]:
spos = bto_tetra.get_scaled_positions()
symbols = bto_tetra.symbols
for i in range(len(spos)):
print(i, symbols[i], sg.equivalent_sites(spos[i]))
0 Ba (array([[0.5 , 0.5 , 0.004136]]), [0])
1 Ti (array([[0. , 0. , 0.524313]]), [0])
2 O (array([[0. , 0.5 , 0.479349],
[0.5 , 0. , 0.479349]]), [0, 0])
3 O (array([[0.5 , 0. , 0.479349],
[0. , 0.5 , 0.479349]]), [0, 0])
4 O (array([[0. , 0. , 0.958854]]), [0])
Finally, we will try to determine the space group for each of the four structures after structural optimization in this tutorial.
The symprec argument of the get_spacegroup method allows you to set a tolerance for the shift of each atom’s coordinates when making the symmetry determination. In this case, the default value of symprec=1e-5 resulted in P 4 m m, so we set a strictly smaller value to determine the space group.
[18]:
symprec = 1e-10
print("atoms1: ", get_spacegroup(atoms1, symprec=symprec).symbol)
print("atoms2: ", get_spacegroup(atoms2, symprec=symprec).symbol)
print("atoms3: ", get_spacegroup(atoms3, symprec=symprec).symbol)
print("atoms4: ", get_spacegroup(atoms4, symprec=symprec).symbol)
atoms1: P 1
atoms2: P 4 m m
atoms3: P 1
atoms4: P 4 m m
Thus, atoms2 and atoms4 with FixSymmetry are P 4 m m preserving the space group, but atoms1 and atoms3 have broken the space group and are now at P 1, which is the least symmetric.
Bravais lattice¶
Now that we have explained the space group, we will explain a related concept, the Bravais lattice.
The space group is a symmetry that depends on both the cell and atomic configuration of the atoms, whereas the Bravais lattice is determined from the cells only. Bravais lattices are known to be one of the following 14 types. These 14 types are also illustrated in the Wikipedia link above.
[19]:
from ase.lattice import CUB, FCC, BCC, TET, BCT, HEX, RHL, ORC, ORCF, ORCI, ORCC, MCL, MCLC, TRI
print(" 1. CUB:", CUB.longname, ",", CUB.crystal_family)
print(" 2. FCC:", FCC.longname, ",", FCC.crystal_family)
print(" 3. BCC:", BCC.longname, ",", BCC.crystal_family)
print(" 4. TET:", TET.longname, ",", TET.crystal_family)
print(" 5. BCT:", BCT.longname, ",", BCT.crystal_family)
print(" 6. HEX:", HEX.longname, ",", HEX.crystal_family)
print(" 7. RHL:", RHL.longname, ",", RHL.crystal_family)
print(" 8. ORC:", ORC.longname, ",", ORC.crystal_family)
print(" 9. ORCF:", ORCF.longname, ",", ORCF.crystal_family)
print("10. ORCI:", ORCI.longname, ",", ORCI.crystal_family)
print("11. ORCC:", ORCC.longname, ",", ORCC.crystal_family)
print("12. MCL :", MCL.longname, ",", MCL.crystal_family)
print("13. MCLC:", MCLC.longname, ",", MCLC.crystal_family)
print("14. TRI :", TRI.longname, ",", TRI.crystal_family)
1. CUB: primitive cubic , cubic
2. FCC: face-centred cubic , cubic
3. BCC: body-centred cubic , cubic
4. TET: primitive tetragonal , tetragonal
5. BCT: body-centred tetragonal , tetragonal
6. HEX: primitive hexagonal , hexagonal
7. RHL: primitive rhombohedral , hexagonal
8. ORC: primitive orthorhombic , orthorhombic
9. ORCF: face-centred orthorhombic , orthorhombic
10. ORCI: body-centred orthorhombic , orthorhombic
11. ORCC: base-centred orthorhombic , orthorhombic
12. MCL : primitive monoclinic , monoclinic
13. MCLC: base-centred monoclinic , monoclinic
14. TRI : primitive triclinic , triclinic
In ASE, the get_bravais_lattice method can be used to determine the bravais lattice.
[20]:
bto_tetra.cell.get_bravais_lattice()
[20]:
TET(a=4.0044570000000003773, c=4.2006360000000002586)
For the tetragonal BaTiO3 structure, TET was obtained as expected.
It is important to note that the cell must be a primitive cell.
When you create a Fe BCC structure and use get_bravais_lattice, if you specify a unit cell, it returns BCC as expected, but if we made a structure with cubic unit cell with cubic=True, CUB is returned.
[21]:
from ase.build import bulk
from pfcc_extras.visualize.ase import view_ase_atoms
fe_bcc = bulk("Fe")
view_ase_atoms(fe_bcc, figsize=(4, 4), title="Fe BCC atoms", scale=100, rotation="0x,0y,0z")
[22]:
# view_ngl(fe_bcc)
[23]:
print("fe_bcc Bravais lattice:", fe_bcc.cell.get_bravais_lattice())
print("fe_bcc cell :", fe_bcc.cell)
fe_bcc Bravais lattice: BCC(a=2.87)
fe_bcc cell : Cell([[-1.435, 1.435, 1.435], [1.435, -1.435, 1.435], [1.435, 1.435, -1.435]])
[24]:
fe_cubic = bulk("Fe", cubic=True)
view_ase_atoms(fe_cubic, figsize=(4, 4), title="Fe BCC cubic=True", scale=100, rotation="0x,0y,0z")
[25]:
# view_ngl(fe_cubic)
[26]:
print("fe_cubic Bravais lattice:", fe_cubic.cell.get_bravais_lattice())
print("fe_cubic cell :", fe_cubic.cell)
fe_cubic Bravais lattice: CUB(a=2.87)
fe_cubic cell : Cell([2.87, 2.87, 2.87])
Reference¶
If you would like to learn more about space groups, the following reference will be helpful.
“Space Groups for Solid State Scientists”, Michael Glazer and Gerald Burns.