Coverage for /builds/debichem-team/python-ase/ase/build/bulk.py: 90.20%

1"""Build crystalline systems""" 

2from math import sqrt 

3from typing import Any 

4 

5from ase.atoms import Atoms 

6from ase.data import atomic_numbers, chemical_symbols, reference_states 

7from ase.symbols import string2symbols 

8from ase.utils import plural 

9 

10 

11def incompatible_cell(*, want, have): 

12 return RuntimeError(f'Cannot create {want} cell for {have} structure') 

13 

14 

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. 

30 

31 Crystal structure and lattice constant(s) will be guessed if not 

32 provided. 

33 

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 """ 

59 

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 

63 

64 if covera is not None and c is not None: 

65 raise ValueError("Don't specify both c and c/a!") 

66 

67 xref = '' 

68 ref: Any = {} 

69 

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'] 

77 

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') 

89 

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} 

101 

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}') 

107 

108 magmom_per_atom = None 

109 if crystalstructure == xref: 

110 magmom_per_atom = ref.get('magmom_per_atom') 

111 

112 if crystalstructure not in structures: 

113 raise ValueError(f'Unknown structure: {crystalstructure}.') 

114 

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)) 

122 

123 if alpha is None: 

124 alpha = ref.get('alpha') 

125 

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}"') 

133 

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 

138 

139 if crystalstructure in ['hcp', 'wurtzite']: 

140 if c is not None: 

141 covera = c / a 

142 elif covera is None: 

143 if xref == crystalstructure: 

144 covera = ref['c/a'] 

145 else: 

146 covera = sqrt(8 / 3) 

147 

148 if covera is None: 

149 covera = ref.get('c/a') 

150 if c is None and covera is not None: 

151 c = covera * a 

152 

153 if crystalstructure == 'bct': 

154 from ase.lattice import BCT 

155 if basis is None: 

156 basis = ref.get('basis') 

157 if basis is not None: 

158 natoms = len(basis) 

159 lat = BCT(a=a, c=c) 

160 atoms = Atoms([name] * natoms, cell=lat.tocell(), pbc=True, 

161 scaled_positions=basis) 

162 elif crystalstructure == 'rhombohedral': 

163 atoms = _build_rhl(name, a, alpha, basis) 

164 elif crystalstructure == 'orthorhombic': 

165 atoms = Atoms(name, cell=[a, b, c], pbc=True) 

166 elif orthorhombic: 

167 atoms = _orthorhombic_bulk(name, crystalstructure, a, covera, u) 

168 elif cubic: 

169 atoms = _cubic_bulk(name, crystalstructure, a) 

170 else: 

171 atoms = _primitive_bulk(name, crystalstructure, a, covera, u) 

172 

173 if magmom_per_atom is not None: 

174 magmoms = [magmom_per_atom] * len(atoms) 

175 atoms.set_initial_magnetic_moments(magmoms) 

176 

177 if cubic or orthorhombic: 

178 assert atoms.cell.orthorhombic 

179 

180 return atoms 

181 

182 

183def _build_rhl(name, a, alpha, basis): 

184 from ase.lattice import RHL 

185 lat = RHL(a, alpha) 

186 cell = lat.tocell() 

187 if basis is None: 

188 # RHL: Given by A&M as scaled coordinates "x" of cell.sum(0): 

189 basis_x = reference_states[atomic_numbers[name]]['basis_x'] 

190 basis = basis_x[:, None].repeat(3, axis=1) 

191 natoms = len(basis) 

192 return Atoms([name] * natoms, cell=cell, scaled_positions=basis, pbc=True) 

193 

194 

195def _orthorhombic_bulk(name, crystalstructure, a, covera=None, u=None): 

196 if crystalstructure in ('sc', 'bcc', 'cesiumchloride'): 

197 atoms = _cubic_bulk(name, crystalstructure, a) 

198 elif crystalstructure == 'fcc': 

199 b = a / sqrt(2) 

200 cell = (b, b, a) 

201 scaled_positions = ((0.0, 0.0, 0.0), (0.5, 0.5, 0.5)) 

202 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions) 

203 elif crystalstructure == 'hcp': 

204 cell = (a, a * sqrt(3), covera * a) 

205 scaled_positions = [ 

206 (0.0, 0 / 6, 0.0), 

207 (0.5, 3 / 6, 0.0), 

208 (0.5, 1 / 6, 0.5), 

209 (0.0, 4 / 6, 0.5), 

210 ] 

211 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions) 

212 elif crystalstructure == 'diamond': 

213 b = a / sqrt(2) 

214 cell = (b, b, a) 

215 scaled_positions = [ 

216 (0.0, 0.0, 0.0), (0.5, 0.0, 0.25), 

217 (0.5, 0.5, 0.5), (0.0, 0.5, 0.75), 

218 ] 

219 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions) 

220 elif crystalstructure == 'rocksalt': 

221 b = a / sqrt(2) 

222 cell = (b, b, a) 

223 scaled_positions = [ 

224 (0.0, 0.0, 0.0), (0.5, 0.5, 0.0), 

225 (0.5, 0.5, 0.5), (0.0, 0.0, 0.5), 

226 ] 

227 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions) 

228 elif crystalstructure == 'zincblende': 

229 symbol0, symbol1 = string2symbols(name) 

230 atoms = _orthorhombic_bulk(symbol0, 'diamond', a) 

231 atoms.symbols[[1, 3]] = symbol1 

232 elif crystalstructure == 'wurtzite': 

233 cell = (a, a * sqrt(3), covera * a) 

234 u = u or 0.25 + 1 / 3 / covera**2 

235 scaled_positions = [ 

236 (0.0, 0 / 6, 0.0), (0.0, 2 / 6, 0.5 - u), 

237 (0.0, 2 / 6, 0.5), (0.0, 0 / 6, 1.0 - u), 

238 (0.5, 3 / 6, 0.0), (0.5, 5 / 6, 0.5 - u), 

239 (0.5, 5 / 6, 0.5), (0.5, 3 / 6, 1.0 - u), 

240 ] 

241 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions) 

242 else: 

243 raise incompatible_cell(want='orthorhombic', have=crystalstructure) 

244 

245 atoms.pbc = True 

246 

247 return atoms 

248 

249 

250def _cubic_bulk(name: str, crystalstructure: str, a: float) -> Atoms: 

251 cell = (a, a, a) 

252 if crystalstructure == 'sc': 

253 atoms = Atoms(name, cell=cell) 

254 elif crystalstructure == 'fcc': 

255 scaled_positions = [ 

256 (0.0, 0.0, 0.0), 

257 (0.0, 0.5, 0.5), 

258 (0.5, 0.0, 0.5), 

259 (0.5, 0.5, 0.0), 

260 ] 

261 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions) 

262 elif crystalstructure == 'bcc': 

263 scaled_positions = [ 

264 (0.0, 0.0, 0.0), 

265 (0.5, 0.5, 0.5), 

266 ] 

267 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions) 

268 elif crystalstructure == 'diamond': 

269 scaled_positions = [ 

270 (0.0, 0.0, 0.0), (0.25, 0.25, 0.25), 

271 (0.0, 0.5, 0.5), (0.25, 0.75, 0.75), 

272 (0.5, 0.0, 0.5), (0.75, 0.25, 0.75), 

273 (0.5, 0.5, 0.0), (0.75, 0.75, 0.25), 

274 ] 

275 atoms = Atoms(8 * name, cell=cell, scaled_positions=scaled_positions) 

276 elif crystalstructure == 'cesiumchloride': 

277 symbol0, symbol1 = string2symbols(name) 

278 atoms = _cubic_bulk(symbol0, 'bcc', a) 

279 atoms.symbols[[1]] = symbol1 

280 elif crystalstructure == 'zincblende': 

281 symbol0, symbol1 = string2symbols(name) 

282 atoms = _cubic_bulk(symbol0, 'diamond', a) 

283 atoms.symbols[[1, 3, 5, 7]] = symbol1 

284 elif crystalstructure == 'rocksalt': 

285 scaled_positions = [ 

286 (0.0, 0.0, 0.0), (0.5, 0.0, 0.0), 

287 (0.0, 0.5, 0.5), (0.5, 0.5, 0.5), 

288 (0.5, 0.0, 0.5), (0.0, 0.0, 0.5), 

289 (0.5, 0.5, 0.0), (0.0, 0.5, 0.0), 

290 ] 

291 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions) 

292 elif crystalstructure == 'fluorite': 

293 scaled_positions = [ 

294 (0.00, 0.00, 0.00), (0.25, 0.25, 0.25), (0.75, 0.75, 0.75), 

295 (0.00, 0.50, 0.50), (0.25, 0.75, 0.75), (0.75, 0.25, 0.25), 

296 (0.50, 0.00, 0.50), (0.75, 0.25, 0.75), (0.25, 0.75, 0.25), 

297 (0.50, 0.50, 0.00), (0.75, 0.75, 0.25), (0.25, 0.25, 0.75), 

298 ] 

299 atoms = Atoms(4 * name, cell=cell, scaled_positions=scaled_positions) 

300 else: 

301 raise incompatible_cell(want='cubic', have=crystalstructure) 

302 

303 atoms.pbc = True 

304 

305 return atoms 

306 

307 

308def _primitive_bulk(name, crystalstructure, a, covera=None, u=None): 

309 if crystalstructure == 'sc': 

310 atoms = Atoms(name, cell=(a, a, a)) 

311 elif crystalstructure == 'fcc': 

312 b = 0.5 * a 

313 cell = ((0, b, b), (b, 0, b), (b, b, 0)) 

314 atoms = Atoms(name, cell=cell) 

315 elif crystalstructure == 'bcc': 

316 b = 0.5 * a 

317 cell = ((-b, b, b), (b, -b, b), (b, b, -b)) 

318 atoms = Atoms(name, cell=cell) 

319 elif crystalstructure == 'hcp': 

320 c = covera * a 

321 cell = ((a, 0, 0), (-0.5 * a, 0.5 * sqrt(3) * a, 0), (0, 0, c)) 

322 scaled_positions = [ 

323 (0 / 3, 0 / 3, 0.0), 

324 (1 / 3, 2 / 3, 0.5), 

325 ] 

326 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions) 

327 elif crystalstructure == 'diamond': 

328 atoms = \ 

329 _primitive_bulk(name, 'fcc', a) + \ 

330 _primitive_bulk(name, 'fcc', a) 

331 atoms.positions[1, :] += 0.25 * a 

332 elif crystalstructure == 'rocksalt': 

333 symbol0, symbol1 = string2symbols(name) 

334 atoms = \ 

335 _primitive_bulk(symbol0, 'fcc', a) + \ 

336 _primitive_bulk(symbol1, 'fcc', a) 

337 atoms.positions[1, 0] += 0.5 * a 

338 elif crystalstructure == 'cesiumchloride': 

339 symbol0, symbol1 = string2symbols(name) 

340 atoms = \ 

341 _primitive_bulk(symbol0, 'sc', a) + \ 

342 _primitive_bulk(symbol1, 'sc', a) 

343 atoms.positions[1, :] += 0.5 * a 

344 elif crystalstructure == 'zincblende': 

345 symbol0, symbol1 = string2symbols(name) 

346 atoms = \ 

347 _primitive_bulk(symbol0, 'fcc', a) + \ 

348 _primitive_bulk(symbol1, 'fcc', a) 

349 atoms.positions[1, :] += 0.25 * a 

350 elif crystalstructure == 'fluorite': 

351 symbol0, symbol1, symbol2 = string2symbols(name) 

352 atoms = \ 

353 _primitive_bulk(symbol0, 'fcc', a) + \ 

354 _primitive_bulk(symbol1, 'fcc', a) + \ 

355 _primitive_bulk(symbol2, 'fcc', a) 

356 atoms.positions[1, :] += 0.25 * a 

357 atoms.positions[2, :] += 0.75 * a 

358 elif crystalstructure == 'wurtzite': 

359 c = covera * a 

360 cell = ((a, 0, 0), (-0.5 * a, 0.5 * sqrt(3) * a, 0), (0, 0, c)) 

361 u = u or 0.25 + 1 / 3 / covera**2 

362 scaled_positions = [ 

363 (0 / 3, 0 / 3, 0.0), (1 / 3, 2 / 3, 0.5 - u), 

364 (1 / 3, 2 / 3, 0.5), (0 / 3, 0 / 3, 1.0 - u), 

365 ] 

366 atoms = Atoms(2 * name, cell=cell, scaled_positions=scaled_positions) 

367 else: 

368 raise incompatible_cell(want='primitive', have=crystalstructure) 

369 

370 atoms.pbc = True 

371 

372 return atoms