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

332 

333 return atoms