Coverage for /builds/ase/ase/ase/build/bulk.py: 90.20%
204 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3"""Build crystalline systems"""
4from math import sqrt
5from typing import Any
7from ase.atoms import Atoms
8from ase.data import atomic_numbers, chemical_symbols, reference_states
9from ase.symbols import string2symbols
10from ase.utils import plural
13def incompatible_cell(*, want, have):
14 return RuntimeError(f'Cannot create {want} cell for {have} structure')
17def bulk(
18 name: str,
19 crystalstructure: str = None,
20 a: float = None,
21 b: float = None,
22 c: float = None,
23 *,
24 alpha: float = None,
25 covera: float = None,
26 u: float = None,
27 orthorhombic: bool = False,
28 cubic: bool = False,
29 basis=None,
30) -> Atoms:
31 """Creating bulk systems.
33 Crystal structure and lattice constant(s) will be guessed if not
34 provided.
36 name: str
37 Chemical symbol or symbols as in 'MgO' or 'NaCl'.
38 crystalstructure: str
39 Must be one of sc, fcc, bcc, tetragonal, bct, hcp, rhombohedral,
40 orthorhombic, mcl, diamond, zincblende, rocksalt, cesiumchloride,
41 fluorite or wurtzite.
42 a: float
43 Lattice constant.
44 b: float
45 Lattice constant. If only a and b is given, b will be interpreted
46 as c instead.
47 c: float
48 Lattice constant.
49 alpha: float
50 Angle in degrees for rhombohedral lattice.
51 covera: float
52 c/a ratio used for hcp. Default is ideal ratio: sqrt(8/3).
53 u: float
54 Internal coordinate for Wurtzite structure.
55 orthorhombic: bool
56 Construct orthorhombic unit cell instead of primitive cell
57 which is the default.
58 cubic: bool
59 Construct cubic unit cell if possible.
60 """
62 if c is None and b is not None:
63 # If user passes (a, b) positionally, we want it as (a, c) instead:
64 c, b = b, c
66 if covera is not None and c is not None:
67 raise ValueError("Don't specify both c and c/a!")
69 xref = ''
70 ref: Any = {}
72 if name in chemical_symbols: # single element
73 atomic_number = atomic_numbers[name]
74 ref = reference_states[atomic_number]
75 if ref is None:
76 ref = {} # easier to 'get' things from empty dictionary than None
77 else:
78 xref = ref['symmetry']
80 if crystalstructure is None:
81 # `ref` requires `basis` but not given and not pre-defined
82 if basis is None and 'basis' in ref and ref['basis'] is None:
83 raise ValueError('This structure requires an atomic basis')
84 if xref == 'cubic':
85 # P and Mn are listed as 'cubic' but the lattice constants
86 # are 7 and 9. They must be something other than simple cubic
87 # then. We used to just return the cubic one but that must
88 # have been wrong somehow. --askhl
89 raise ValueError(
90 f'The reference structure of {name} is not implemented')
92 # Mapping of name to number of atoms in primitive cell.
93 structures = {'sc': 1, 'fcc': 1, 'bcc': 1,
94 'tetragonal': 1,
95 'bct': 1,
96 'hcp': 1,
97 'rhombohedral': 1,
98 'orthorhombic': 1,
99 'mcl': 1,
100 'diamond': 1,
101 'zincblende': 2, 'rocksalt': 2, 'cesiumchloride': 2,
102 'fluorite': 3, 'wurtzite': 2}
104 if crystalstructure is None:
105 crystalstructure = xref
106 if crystalstructure not in structures:
107 raise ValueError(f'No suitable reference data for bulk {name}.'
108 f' Reference data: {ref}')
110 magmom_per_atom = None
111 if crystalstructure == xref:
112 magmom_per_atom = ref.get('magmom_per_atom')
114 if crystalstructure not in structures:
115 raise ValueError(f'Unknown structure: {crystalstructure}.')
117 # Check name:
118 natoms = len(string2symbols(name))
119 natoms0 = structures[crystalstructure]
120 if natoms != natoms0:
121 raise ValueError('Please specify {} for {} and not {}'
122 .format(plural(natoms0, 'atom'),
123 crystalstructure, natoms))
125 if alpha is None:
126 alpha = ref.get('alpha')
128 if a is None:
129 if xref != crystalstructure:
130 raise ValueError('You need to specify the lattice constant.')
131 if 'a' in ref:
132 a = ref['a']
133 else:
134 raise KeyError(f'No reference lattice parameter "a" for "{name}"')
136 if b is None:
137 bovera = ref.get('b/a')
138 if bovera is not None and a is not None:
139 b = bovera * a
141 if crystalstructure in ['hcp', 'wurtzite']:
142 if c is not None:
143 covera = c / a
144 elif covera is None:
145 if xref == crystalstructure:
146 covera = ref['c/a']
147 else:
148 covera = sqrt(8 / 3)
150 if covera is None:
151 covera = ref.get('c/a')
152 if c is None and covera is not None:
153 c = covera * a
155 if crystalstructure == 'bct':
156 from ase.lattice import BCT
157 if basis is None:
158 basis = ref.get('basis')
159 if basis is not None:
160 natoms = len(basis)
161 lat = BCT(a=a, c=c)
162 atoms = Atoms([name] * natoms, cell=lat.tocell(), pbc=True,
163 scaled_positions=basis)
164 elif crystalstructure == 'rhombohedral':
165 atoms = _build_rhl(name, a, alpha, basis)
166 elif crystalstructure == 'orthorhombic':
167 atoms = Atoms(name, cell=[a, b, c], pbc=True)
168 elif orthorhombic:
169 atoms = _orthorhombic_bulk(name, crystalstructure, a, covera, u)
170 elif cubic:
171 atoms = _cubic_bulk(name, crystalstructure, a)
172 else:
173 atoms = _primitive_bulk(name, crystalstructure, a, covera, u)
175 if magmom_per_atom is not None:
176 magmoms = [magmom_per_atom] * len(atoms)
177 atoms.set_initial_magnetic_moments(magmoms)
179 if cubic or orthorhombic:
180 assert atoms.cell.orthorhombic
182 return atoms
185def _build_rhl(name, a, alpha, basis):
186 from ase.lattice import RHL
187 lat = RHL(a, alpha)
188 cell = lat.tocell()
189 if basis is None:
190 # RHL: Given by A&M as scaled coordinates "x" of cell.sum(0):
191 basis_x = reference_states[atomic_numbers[name]]['basis_x']
192 basis = basis_x[:, None].repeat(3, axis=1)
193 natoms = len(basis)
194 return Atoms([name] * natoms, cell=cell, scaled_positions=basis, pbc=True)
197def _orthorhombic_bulk(name, crystalstructure, a, covera=None, u=None):
198 if crystalstructure in ('sc', 'bcc', 'cesiumchloride'):
199 atoms = _cubic_bulk(name, crystalstructure, a)
200 elif crystalstructure == 'fcc':
201 b = a / sqrt(2)
202 cell = (b, b, a)
203 scaled_positions = ((0.0, 0.0, 0.0), (0.5, 0.5, 0.5))
204 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions)
205 elif crystalstructure == 'hcp':
206 cell = (a, a * sqrt(3), covera * a)
207 scaled_positions = [
208 (0.0, 0 / 6, 0.0),
209 (0.5, 3 / 6, 0.0),
210 (0.5, 1 / 6, 0.5),
211 (0.0, 4 / 6, 0.5),
212 ]
213 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions)
214 elif crystalstructure == 'diamond':
215 b = a / sqrt(2)
216 cell = (b, b, a)
217 scaled_positions = [
218 (0.0, 0.0, 0.0), (0.5, 0.0, 0.25),
219 (0.5, 0.5, 0.5), (0.0, 0.5, 0.75),
220 ]
221 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions)
222 elif crystalstructure == 'rocksalt':
223 b = a / sqrt(2)
224 cell = (b, b, a)
225 scaled_positions = [
226 (0.0, 0.0, 0.0), (0.5, 0.5, 0.0),
227 (0.5, 0.5, 0.5), (0.0, 0.0, 0.5),
228 ]
229 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions)
230 elif crystalstructure == 'zincblende':
231 symbol0, symbol1 = string2symbols(name)
232 atoms = _orthorhombic_bulk(symbol0, 'diamond', a)
233 atoms.symbols[[1, 3]] = symbol1
234 elif crystalstructure == 'wurtzite':
235 cell = (a, a * sqrt(3), covera * a)
236 u = u or 0.25 + 1 / 3 / covera**2
237 scaled_positions = [
238 (0.0, 0 / 6, 0.0), (0.0, 2 / 6, 0.5 - u),
239 (0.0, 2 / 6, 0.5), (0.0, 0 / 6, 1.0 - u),
240 (0.5, 3 / 6, 0.0), (0.5, 5 / 6, 0.5 - u),
241 (0.5, 5 / 6, 0.5), (0.5, 3 / 6, 1.0 - u),
242 ]
243 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions)
244 else:
245 raise incompatible_cell(want='orthorhombic', have=crystalstructure)
247 atoms.pbc = True
249 return atoms
252def _cubic_bulk(name: str, crystalstructure: str, a: float) -> Atoms:
253 cell = (a, a, a)
254 if crystalstructure == 'sc':
255 atoms = Atoms(name, cell=cell)
256 elif crystalstructure == 'fcc':
257 scaled_positions = [
258 (0.0, 0.0, 0.0),
259 (0.0, 0.5, 0.5),
260 (0.5, 0.0, 0.5),
261 (0.5, 0.5, 0.0),
262 ]
263 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions)
264 elif crystalstructure == 'bcc':
265 scaled_positions = [
266 (0.0, 0.0, 0.0),
267 (0.5, 0.5, 0.5),
268 ]
269 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions)
270 elif crystalstructure == 'diamond':
271 scaled_positions = [
272 (0.0, 0.0, 0.0), (0.25, 0.25, 0.25),
273 (0.0, 0.5, 0.5), (0.25, 0.75, 0.75),
274 (0.5, 0.0, 0.5), (0.75, 0.25, 0.75),
275 (0.5, 0.5, 0.0), (0.75, 0.75, 0.25),
276 ]
277 atoms = Atoms(8 * name, cell=cell, scaled_positions=scaled_positions)
278 elif crystalstructure == 'cesiumchloride':
279 symbol0, symbol1 = string2symbols(name)
280 atoms = _cubic_bulk(symbol0, 'bcc', a)
281 atoms.symbols[[1]] = symbol1
282 elif crystalstructure == 'zincblende':
283 symbol0, symbol1 = string2symbols(name)
284 atoms = _cubic_bulk(symbol0, 'diamond', a)
285 atoms.symbols[[1, 3, 5, 7]] = symbol1
286 elif crystalstructure == 'rocksalt':
287 scaled_positions = [
288 (0.0, 0.0, 0.0), (0.5, 0.0, 0.0),
289 (0.0, 0.5, 0.5), (0.5, 0.5, 0.5),
290 (0.5, 0.0, 0.5), (0.0, 0.0, 0.5),
291 (0.5, 0.5, 0.0), (0.0, 0.5, 0.0),
292 ]
293 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions)
294 elif crystalstructure == 'fluorite':
295 scaled_positions = [
296 (0.00, 0.00, 0.00), (0.25, 0.25, 0.25), (0.75, 0.75, 0.75),
297 (0.00, 0.50, 0.50), (0.25, 0.75, 0.75), (0.75, 0.25, 0.25),
298 (0.50, 0.00, 0.50), (0.75, 0.25, 0.75), (0.25, 0.75, 0.25),
299 (0.50, 0.50, 0.00), (0.75, 0.75, 0.25), (0.25, 0.25, 0.75),
300 ]
301 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions)
302 else:
303 raise incompatible_cell(want='cubic', have=crystalstructure)
305 atoms.pbc = True
307 return atoms
310def _primitive_bulk(name, crystalstructure, a, covera=None, u=None):
311 if crystalstructure == 'sc':
312 atoms = Atoms(name, cell=(a, a, a))
313 elif crystalstructure == 'fcc':
314 b = 0.5 * a
315 cell = ((0, b, b), (b, 0, b), (b, b, 0))
316 atoms = Atoms(name, cell=cell)
317 elif crystalstructure == 'bcc':
318 b = 0.5 * a
319 cell = ((-b, b, b), (b, -b, b), (b, b, -b))
320 atoms = Atoms(name, cell=cell)
321 elif crystalstructure == 'hcp':
322 c = covera * a
323 cell = ((a, 0, 0), (-0.5 * a, 0.5 * sqrt(3) * a, 0), (0, 0, c))
324 scaled_positions = [
325 (0 / 3, 0 / 3, 0.0),
326 (1 / 3, 2 / 3, 0.5),
327 ]
328 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions)
329 elif crystalstructure == 'diamond':
330 atoms = \
331 _primitive_bulk(name, 'fcc', a) + \
332 _primitive_bulk(name, 'fcc', a)
333 atoms.positions[1, :] += 0.25 * a
334 elif crystalstructure == 'rocksalt':
335 symbol0, symbol1 = string2symbols(name)
336 atoms = \
337 _primitive_bulk(symbol0, 'fcc', a) + \
338 _primitive_bulk(symbol1, 'fcc', a)
339 atoms.positions[1, 0] += 0.5 * a
340 elif crystalstructure == 'cesiumchloride':
341 symbol0, symbol1 = string2symbols(name)
342 atoms = \
343 _primitive_bulk(symbol0, 'sc', a) + \
344 _primitive_bulk(symbol1, 'sc', a)
345 atoms.positions[1, :] += 0.5 * a
346 elif crystalstructure == 'zincblende':
347 symbol0, symbol1 = string2symbols(name)
348 atoms = \
349 _primitive_bulk(symbol0, 'fcc', a) + \
350 _primitive_bulk(symbol1, 'fcc', a)
351 atoms.positions[1, :] += 0.25 * a
352 elif crystalstructure == 'fluorite':
353 symbol0, symbol1, symbol2 = string2symbols(name)
354 atoms = \
355 _primitive_bulk(symbol0, 'fcc', a) + \
356 _primitive_bulk(symbol1, 'fcc', a) + \
357 _primitive_bulk(symbol2, 'fcc', a)
358 atoms.positions[1, :] += 0.25 * a
359 atoms.positions[2, :] += 0.75 * a
360 elif crystalstructure == 'wurtzite':
361 c = covera * a
362 cell = ((a, 0, 0), (-0.5 * a, 0.5 * sqrt(3) * a, 0), (0, 0, c))
363 u = u or 0.25 + 1 / 3 / covera**2
364 scaled_positions = [
365 (0 / 3, 0 / 3, 0.0), (1 / 3, 2 / 3, 0.5 - u),
366 (1 / 3, 2 / 3, 0.5), (0 / 3, 0 / 3, 1.0 - u),
367 ]
368 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions)
369 else:
370 raise incompatible_cell(want='primitive', have=crystalstructure)
372 atoms.pbc = True
374 return atoms