Coverage for /builds/ase/ase/ase/build/bulk.py : 87.72%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1from math import sqrt
3from ase.atoms import Atoms
4from ase.symbols import string2symbols
5from ase.data import reference_states, atomic_numbers, chemical_symbols
6from ase.utils import plural
9def incompatible_cell(*, want, have):
10 return RuntimeError('Cannot create {} cell for {} structure'
11 .format(want, have))
14def bulk(name, crystalstructure=None, a=None, b=None, c=None, *, alpha=None,
15 covera=None, u=None, orthorhombic=False, cubic=False,
16 basis=None):
17 """Creating bulk systems.
19 Crystal structure and lattice constant(s) will be guessed if not
20 provided.
22 name: str
23 Chemical symbol or symbols as in 'MgO' or 'NaCl'.
24 crystalstructure: str
25 Must be one of sc, fcc, bcc, tetragonal, bct, hcp, rhombohedral,
26 orthorhombic, mcl, diamond, zincblende, rocksalt, cesiumchloride,
27 fluorite or wurtzite.
28 a: float
29 Lattice constant.
30 b: float
31 Lattice constant. If only a and b is given, b will be interpreted
32 as c instead.
33 c: float
34 Lattice constant.
35 alpha: float
36 Angle in degrees for rhombohedral lattice.
37 covera: float
38 c/a ratio used for hcp. Default is ideal ratio: sqrt(8/3).
39 u: float
40 Internal coordinate for Wurtzite structure.
41 orthorhombic: bool
42 Construct orthorhombic unit cell instead of primitive cell
43 which is the default.
44 cubic: bool
45 Construct cubic unit cell if possible.
46 """
48 if c is None and b is not None:
49 # If user passes (a, b) positionally, we want it as (a, c) instead:
50 c, b = b, c
52 if covera is not None and c is not None:
53 raise ValueError("Don't specify both c and c/a!")
55 xref = None
56 ref = {}
58 if name in chemical_symbols:
59 Z = atomic_numbers[name]
60 ref = reference_states[Z]
61 if ref is not None:
62 xref = ref['symmetry']
64 # If user did not specify crystal structure, and no basis
65 # is given, and the reference state says we need one, but
66 # does not have one, then we can't proceed.
67 if (crystalstructure is None and basis is None
68 and 'basis' in ref and ref['basis'] is None):
69 # XXX This is getting much too complicated, we need to split
70 # this function up. A lot.
71 raise RuntimeError('This structure requires an atomic basis')
73 if ref is None:
74 ref = {} # easier to 'get' things from empty dictionary than None
76 if xref == 'cubic':
77 # P and Mn are listed as 'cubic' but the lattice constants
78 # are 7 and 9. They must be something other than simple cubic
79 # then. We used to just return the cubic one but that must
80 # have been wrong somehow. --askhl
81 raise RuntimeError('Only simple cubic ("sc") supported')
83 # Mapping of name to number of atoms in primitive cell.
84 structures = {'sc': 1, 'fcc': 1, 'bcc': 1,
85 'tetragonal': 1,
86 'bct': 1,
87 'hcp': 1,
88 'rhombohedral': 1,
89 'orthorhombic': 1,
90 'mcl': 1,
91 'diamond': 1,
92 'zincblende': 2, 'rocksalt': 2, 'cesiumchloride': 2,
93 'fluorite': 3, 'wurtzite': 2}
95 if crystalstructure is None:
96 crystalstructure = xref
97 if crystalstructure not in structures:
98 raise ValueError('No suitable reference data for bulk {}.'
99 ' Reference data: {}'
100 .format(name, ref))
102 magmom_per_atom = None
103 if crystalstructure == xref:
104 magmom_per_atom = ref.get('magmom_per_atom')
106 if crystalstructure not in structures:
107 raise ValueError('Unknown structure: {}.'
108 .format(crystalstructure))
110 # Check name:
111 natoms = len(string2symbols(name))
112 natoms0 = structures[crystalstructure]
113 if natoms != natoms0:
114 raise ValueError('Please specify {} for {} and not {}'
115 .format(plural(natoms0, 'atom'),
116 crystalstructure, natoms))
118 if alpha is None:
119 alpha = ref.get('alpha')
121 if a is None:
122 if xref != crystalstructure:
123 raise ValueError('You need to specify the lattice constant.')
124 try:
125 a = ref['a']
126 except KeyError:
127 raise KeyError('No reference lattice parameter "a" for "{}"'
128 .format(name))
130 if b is None:
131 bovera = ref.get('b/a')
132 if bovera is not None and a is not None:
133 b = bovera * a
135 if crystalstructure in ['hcp', 'wurtzite']:
136 if cubic:
137 raise incompatible_cell(want='cubic', have=crystalstructure)
139 if c is not None:
140 covera = c / a
141 elif covera is None:
142 if xref == crystalstructure:
143 covera = ref['c/a']
144 else:
145 covera = sqrt(8 / 3)
147 if covera is None:
148 covera = ref.get('c/a')
149 if c is None and covera is not None:
150 c = covera * a
152 if orthorhombic and crystalstructure not in ['sc', 'tetragonal',
153 'orthorhombic']:
154 atoms = _orthorhombic_bulk(name, crystalstructure, a, covera, u)
155 elif cubic and crystalstructure in ['bcc', 'cesiumchloride']:
156 atoms = _orthorhombic_bulk(name, crystalstructure, a, covera)
157 elif cubic and crystalstructure != 'sc':
158 atoms = _cubic_bulk(name, crystalstructure, a)
159 elif crystalstructure == 'sc':
160 atoms = Atoms(name, cell=(a, a, a), pbc=True)
161 elif crystalstructure == 'fcc':
162 b = a / 2
163 atoms = Atoms(name, cell=[(0, b, b), (b, 0, b), (b, b, 0)], pbc=True)
164 elif crystalstructure == 'bcc':
165 b = a / 2
166 atoms = Atoms(name, cell=[(-b, b, b), (b, -b, b), (b, b, -b)],
167 pbc=True)
168 elif crystalstructure == 'hcp':
169 atoms = Atoms(2 * name,
170 scaled_positions=[(0, 0, 0),
171 (1 / 3, 2 / 3, 0.5)],
172 cell=[(a, 0, 0),
173 (-a / 2, a * sqrt(3) / 2, 0),
174 (0, 0, covera * a)],
175 pbc=True)
176 elif crystalstructure == 'diamond':
177 atoms = bulk(2 * name, 'zincblende', a)
178 elif crystalstructure == 'zincblende':
179 s1, s2 = string2symbols(name)
180 atoms = bulk(s1, 'fcc', a) + bulk(s2, 'fcc', a)
181 atoms.positions[1] += a / 4
182 elif crystalstructure == 'rocksalt':
183 s1, s2 = string2symbols(name)
184 atoms = bulk(s1, 'fcc', a) + bulk(s2, 'fcc', a)
185 atoms.positions[1, 0] += a / 2
186 elif crystalstructure == 'cesiumchloride':
187 s1, s2 = string2symbols(name)
188 atoms = bulk(s1, 'sc', a) + bulk(s2, 'sc', a)
189 atoms.positions[1, :] += a / 2
190 elif crystalstructure == 'fluorite':
191 s1, s2, s3 = string2symbols(name)
192 atoms = bulk(s1, 'fcc', a) + bulk(s2, 'fcc', a) + bulk(s3, 'fcc', a)
193 atoms.positions[1, :] += a / 4
194 atoms.positions[2, :] += a * 3 / 4
195 elif crystalstructure == 'wurtzite':
196 u = u or 0.25 + 1 / 3 / covera**2
197 atoms = Atoms(2 * name,
198 scaled_positions=[(0, 0, 0),
199 (1 / 3, 2 / 3, 0.5 - u),
200 (1 / 3, 2 / 3, 0.5),
201 (0, 0, 1 - u)],
202 cell=[(a, 0, 0),
203 (-a / 2, a * sqrt(3) / 2, 0),
204 (0, 0, a * covera)],
205 pbc=True)
206 elif crystalstructure == 'bct':
207 from ase.lattice import BCT
208 if basis is None:
209 basis = ref.get('basis')
210 if basis is not None:
211 natoms = len(basis)
212 lat = BCT(a=a, c=c)
213 atoms = Atoms([name] * natoms, cell=lat.tocell(), pbc=True,
214 scaled_positions=basis)
215 elif crystalstructure == 'rhombohedral':
216 atoms = _build_rhl(name, a, alpha, basis)
217 elif crystalstructure == 'orthorhombic':
218 atoms = Atoms(name, cell=[a, b, c], pbc=True)
219 else:
220 raise ValueError(f'Unknown crystal structure: {crystalstructure!r}')
222 if magmom_per_atom is not None:
223 magmoms = [magmom_per_atom] * len(atoms)
224 atoms.set_initial_magnetic_moments(magmoms)
226 if orthorhombic:
227 assert atoms.cell.orthorhombic
229 if cubic:
230 assert abs(atoms.cell.angles() - 90).all() < 1e-10
232 return atoms
235def _build_rhl(name, a, alpha, basis):
236 from ase.lattice import RHL
237 lat = RHL(a, alpha)
238 cell = lat.tocell()
239 if basis is None:
240 # RHL: Given by A&M as scaled coordinates "x" of cell.sum(0):
241 basis_x = reference_states[atomic_numbers[name]]['basis_x']
242 basis = basis_x[:, None].repeat(3, axis=1)
243 natoms = len(basis)
244 return Atoms([name] * natoms, cell=cell, scaled_positions=basis, pbc=True)
247def _orthorhombic_bulk(name, crystalstructure, a, covera=None, u=None):
248 if crystalstructure == 'fcc':
249 b = a / sqrt(2)
250 atoms = Atoms(2 * name, cell=(b, b, a), pbc=True,
251 scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)])
252 elif crystalstructure == 'bcc':
253 atoms = Atoms(2 * name, cell=(a, a, a), pbc=True,
254 scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)])
255 elif crystalstructure == 'hcp':
256 atoms = Atoms(4 * name,
257 cell=(a, a * sqrt(3), covera * a),
258 scaled_positions=[(0, 0, 0),
259 (0.5, 0.5, 0),
260 (0.5, 1 / 6, 0.5),
261 (0, 2 / 3, 0.5)],
262 pbc=True)
263 elif crystalstructure == 'diamond':
264 atoms = _orthorhombic_bulk(2 * name, 'zincblende', a)
265 elif crystalstructure == 'zincblende':
266 s1, s2 = string2symbols(name)
267 b = a / sqrt(2)
268 atoms = Atoms(2 * name, cell=(b, b, a), pbc=True,
269 scaled_positions=[(0, 0, 0), (0.5, 0, 0.25),
270 (0.5, 0.5, 0.5), (0, 0.5, 0.75)])
271 elif crystalstructure == 'rocksalt':
272 s1, s2 = string2symbols(name)
273 b = a / sqrt(2)
274 atoms = Atoms(2 * name, cell=(b, b, a), pbc=True,
275 scaled_positions=[(0, 0, 0), (0.5, 0.5, 0),
276 (0.5, 0.5, 0.5), (0, 0, 0.5)])
277 elif crystalstructure == 'cesiumchloride':
278 atoms = Atoms(name, cell=(a, a, a), pbc=True,
279 scaled_positions=[(0, 0, 0), (0.5, 0.5, 0.5)])
280 elif crystalstructure == 'wurtzite':
281 u = u or 0.25 + 1 / 3 / covera**2
282 atoms = Atoms(4 * name,
283 cell=(a, a * 3**0.5, covera * a),
284 scaled_positions=[(0, 0, 0),
285 (0, 1 / 3, 0.5 - u),
286 (0, 1 / 3, 0.5),
287 (0, 0, 1 - u),
288 (0.5, 0.5, 0),
289 (0.5, 5 / 6, 0.5 - u),
290 (0.5, 5 / 6, 0.5),
291 (0.5, 0.5, 1 - u)],
292 pbc=True)
293 else:
294 raise incompatible_cell(want='orthorhombic', have=crystalstructure)
296 return atoms
299def _cubic_bulk(name, crystalstructure, a):
300 if crystalstructure == 'fcc':
301 atoms = Atoms(4 * name, cell=(a, a, a), pbc=True,
302 scaled_positions=[(0, 0, 0), (0, 0.5, 0.5),
303 (0.5, 0, 0.5), (0.5, 0.5, 0)])
304 elif crystalstructure == 'diamond':
305 atoms = _cubic_bulk(2 * name, 'zincblende', a)
306 elif crystalstructure == 'zincblende':
307 atoms = Atoms(4 * name, cell=(a, a, a), pbc=True,
308 scaled_positions=[(0, 0, 0), (0.25, 0.25, 0.25),
309 (0, 0.5, 0.5), (0.25, 0.75, 0.75),
310 (0.5, 0, 0.5), (0.75, 0.25, 0.75),
311 (0.5, 0.5, 0), (0.75, 0.75, 0.25)])
312 elif crystalstructure == 'rocksalt':
313 atoms = Atoms(4 * name, cell=(a, a, a), pbc=True,
314 scaled_positions=[(0, 0, 0), (0.5, 0, 0),
315 (0, 0.5, 0.5), (0.5, 0.5, 0.5),
316 (0.5, 0, 0.5), (0, 0, 0.5),
317 (0.5, 0.5, 0), (0, 0.5, 0)])
318 elif crystalstructure == 'fluorite':
319 atoms = Atoms(
320 4 * name,
321 cell=(a, a, a),
322 pbc=True,
323 scaled_positions=[
324 (0.00, 0.00, 0.00), (0.25, 0.25, 0.25), (0.75, 0.75, 0.75),
325 (0.00, 0.50, 0.50), (0.25, 0.75, 0.75), (0.75, 0.25, 0.25),
326 (0.50, 0.00, 0.50), (0.75, 0.25, 0.75), (0.25, 0.75, 0.25),
327 (0.50, 0.50, 0.00), (0.75, 0.75, 0.25), (0.25, 0.25, 0.75),
328 ],
329 )
330 else:
331 raise incompatible_cell(want='cubic', have=crystalstructure)
333 return atoms