Hide keyboard shortcuts

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 

2 

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 

7 

8 

9def incompatible_cell(*, want, have): 

10 return RuntimeError('Cannot create {} cell for {} structure' 

11 .format(want, have)) 

12 

13 

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. 

18 

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

20 provided. 

21 

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

47 

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 

51 

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

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

54 

55 xref = None 

56 ref = {} 

57 

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

63 

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

72 

73 if ref is None: 

74 ref = {} # easier to 'get' things from empty dictionary than None 

75 

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

82 

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} 

94 

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

101 

102 magmom_per_atom = None 

103 if crystalstructure == xref: 

104 magmom_per_atom = ref.get('magmom_per_atom') 

105 

106 if crystalstructure not in structures: 

107 raise ValueError('Unknown structure: {}.' 

108 .format(crystalstructure)) 

109 

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

117 

118 if alpha is None: 

119 alpha = ref.get('alpha') 

120 

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

129 

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 

134 

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

136 if cubic: 

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

138 

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) 

146 

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 

151 

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

221 

222 if magmom_per_atom is not None: 

223 magmoms = [magmom_per_atom] * len(atoms) 

224 atoms.set_initial_magnetic_moments(magmoms) 

225 

226 if orthorhombic: 

227 assert atoms.cell.orthorhombic 

228 

229 if cubic: 

230 assert abs(atoms.cell.angles() - 90).all() < 1e-10 

231 

232 return atoms 

233 

234 

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) 

245 

246 

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) 

295 

296 return atoms 

297 

298 

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) 

320 

321 return atoms