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