Coverage for /builds/ase/ase/ase/build/bulk.py : 82.25%
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 else:
319 raise incompatible_cell(want='cubic', have=crystalstructure)
321 return atoms