Coverage for ase / spacegroup / xtal.py: 93.33%
75 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +0000
1# fmt: off
3# Copyright (C) 2010, Jesper Friis
4# (see accompanying license files for details).
6"""
7A module for ASE for simple creation of crystalline structures from
8knowledge of the space group.
10"""
12from typing import Any
14import numpy as np
15from scipy import spatial
17import ase
18from ase.geometry import cellpar_to_cell
19from ase.spacegroup import Spacegroup
20from ase.symbols import string2symbols
22__all__ = ['crystal']
25def crystal(symbols=None, basis=None, occupancies=None, spacegroup=1, setting=1,
26 cell=None, cellpar=None,
27 ab_normal=(0, 0, 1), a_direction=None, size=(1, 1, 1),
28 onduplicates='warn', symprec=0.001,
29 pbc=True, primitive_cell=False, **kwargs) -> ase.Atoms:
30 """Create an Atoms instance for a conventional unit cell of a
31 space group.
33 Parameters
34 ----------
36 symbols : str | sequence of str | sequence of Atom | Atoms
37 Element symbols of the unique sites. Can either be a string
38 formula or a sequence of element symbols. E.g. ('Na', 'Cl')
39 and 'NaCl' are equivalent. Can also be given as a sequence of
40 Atom objects or an Atoms object.
41 basis : list of scaled coordinates
42 Positions of the unique sites corresponding to symbols given
43 either as scaled positions or through an atoms instance. Not
44 needed if *symbols* is a sequence of Atom objects or an Atoms
45 object.
46 occupancies : list of site occupancies
47 Occupancies of the unique sites. Defaults to 1.0 and thus no mixed
48 occupancies are considered if not explicitly asked for. If occupancies
49 are given, the most dominant species will yield the atomic number.
50 The occupancies in the atoms.info['occupancy'] dictionary will have
51 integers keys converted to strings. The conversion is done in order
52 to avoid unexpected conversions when using the JSON serializer.
53 spacegroup : int | string | Spacegroup instance
54 Space group given either as its number in International Tables
55 or as its Hermann-Mauguin symbol.
56 setting : 1 | 2
57 Space group setting.
58 cell : 3x3 matrix
59 Unit cell vectors.
60 cellpar : [a, b, c, alpha, beta, gamma]
61 Cell parameters with angles in degree. Is not used when `cell`
62 is given.
63 ab_normal : vector
64 Is used to define the orientation of the unit cell relative
65 to the Cartesian system when `cell` is not given. It is the
66 normal vector of the plane spanned by a and b.
67 a_direction : vector
68 Defines the orientation of the unit cell a vector. a will be
69 parallel to the projection of `a_direction` onto the a-b plane.
70 size : 3 positive integers
71 How many times the conventional unit cell should be repeated
72 in each direction.
73 onduplicates : 'keep' | 'replace' | 'warn' | 'error'
74 Action if `basis` contain symmetry-equivalent positions:
75 'keep' - ignore additional symmetry-equivalent positions
76 'replace' - replace
77 'warn' - like 'keep', but issue an UserWarning
78 'error' - raises a SpacegroupValueError
79 symprec : float
80 Minimum "distance" betweed two sites in scaled coordinates
81 before they are counted as the same site.
82 pbc : one or three bools
83 Periodic boundary conditions flags. Examples: True,
84 False, 0, 1, (1, 1, 0), (True, False, False). Default
85 is True.
86 primitive_cell : bool
87 Whether to return the primitive instead of the conventional
88 unit cell.
90 Keyword arguments:
92 All additional keyword arguments are passed on to the Atoms
93 constructor. Currently, probably the most useful additional
94 keyword arguments are `info`, `constraint` and `calculator`.
96 Examples
97 --------
99 Two diamond unit cells (space group number 227)
101 >>> diamond = crystal('C', [(0,0,0)], spacegroup=227,
102 ... cellpar=[3.57, 3.57, 3.57, 90, 90, 90], size=(2,1,1))
103 >>> ase.view(diamond) # doctest: +SKIP
105 A CoSb3 skutterudite unit cell containing 32 atoms
107 >>> skutterudite = crystal(('Co', 'Sb'),
108 ... basis=[(0.25,0.25,0.25), (0.0, 0.335, 0.158)],
109 ... spacegroup=204, cellpar=[9.04, 9.04, 9.04, 90, 90, 90])
110 >>> len(skutterudite)
111 32
112 """
113 sg = Spacegroup(spacegroup, setting)
114 if (
115 not isinstance(symbols, str) and
116 hasattr(symbols, '__getitem__') and
117 len(symbols) > 0 and
118 isinstance(symbols[0], ase.Atom)):
119 symbols = ase.Atoms(symbols)
120 if isinstance(symbols, ase.Atoms):
121 basis = symbols
122 symbols = basis.get_chemical_symbols()
123 if isinstance(basis, ase.Atoms):
124 basis_coords = basis.get_scaled_positions()
125 if cell is None and cellpar is None:
126 cell = basis.cell
127 if symbols is None:
128 symbols = basis.get_chemical_symbols()
129 else:
130 basis_coords = np.array(basis, dtype=float, ndmin=2)
132 if occupancies is not None:
133 occupancies_dict = {}
135 for index, coord in enumerate(basis_coords):
136 # Compute all distances and get indices of nearest atoms
137 dist = spatial.distance.cdist(coord.reshape(1, 3), basis_coords)
138 indices_dist = np.flatnonzero(dist < symprec)
140 occ = {symbols[index]: occupancies[index]}
142 # Check nearest and update occupancy
143 for index_dist in indices_dist:
144 if index == index_dist:
145 continue
146 else:
147 occ.update({symbols[index_dist]: occupancies[index_dist]})
149 occupancies_dict[str(index)] = occ.copy()
151 sites, kinds = sg.equivalent_sites(basis_coords,
152 onduplicates=onduplicates,
153 symprec=symprec)
155 # this is needed to handle deuterium masses
156 masses = None
157 if 'masses' in kwargs:
158 masses = kwargs['masses'][kinds]
159 del kwargs['masses']
161 symbols = parse_symbols(symbols)
163 if occupancies is None:
164 symbols = [symbols[i] for i in kinds]
165 else:
166 # make sure that we put the dominant species there
167 symbols = [sorted(occupancies_dict[str(i)].items(),
168 key=lambda x: x[1])[-1][0] for i in kinds]
170 if cell is None:
171 cell = cellpar_to_cell(cellpar, ab_normal, a_direction)
173 info: dict[str, Any] = {}
174 info['spacegroup'] = sg
175 if primitive_cell:
176 info['unit_cell'] = 'primitive'
177 else:
178 info['unit_cell'] = 'conventional'
180 if 'info' in kwargs:
181 info.update(kwargs['info'])
183 if occupancies is not None:
184 info['occupancy'] = occupancies_dict
186 kwargs['info'] = info
188 atoms = ase.Atoms(symbols,
189 scaled_positions=sites,
190 cell=cell,
191 pbc=pbc,
192 masses=masses,
193 **kwargs)
195 if isinstance(basis, ase.Atoms):
196 for name in basis.arrays:
197 if not atoms.has(name):
198 array = basis.get_array(name)
199 atoms.new_array(name, [array[i] for i in kinds],
200 dtype=array.dtype, shape=array.shape[1:])
202 if kinds:
203 atoms.new_array('spacegroup_kinds', np.asarray(kinds, dtype=int))
205 if primitive_cell:
206 from ase.build import cut
207 prim_cell = sg.scaled_primitive_cell
209 # Preserve calculator if present:
210 calc = atoms.calc
211 atoms = cut(atoms, a=prim_cell[0], b=prim_cell[1], c=prim_cell[2])
212 atoms.calc = calc
214 if size != (1, 1, 1):
215 atoms = atoms.repeat(size)
216 return atoms
219def parse_symbols(symbols):
220 """Return `sumbols` as a sequence of element symbols."""
221 if isinstance(symbols, str):
222 symbols = string2symbols(symbols)
223 return symbols