1# Copyright (C) 2010, Jesper Friis 

2# (see accompanying license files for details). 

3 

4""" 

5A module for ASE for simple creation of crystalline structures from 

6knowledge of the space group. 

7 

8""" 

9 

10from typing import Dict, Any 

11 

12import numpy as np 

13from scipy import spatial 

14 

15import ase 

16from ase.symbols import string2symbols 

17from ase.spacegroup import Spacegroup 

18from ase.geometry import cellpar_to_cell 

19 

20__all__ = ['crystal'] 

21 

22 

23def crystal(symbols=None, basis=None, occupancies=None, spacegroup=1, setting=1, 

24 cell=None, cellpar=None, 

25 ab_normal=(0, 0, 1), a_direction=None, size=(1, 1, 1), 

26 onduplicates='warn', symprec=0.001, 

27 pbc=True, primitive_cell=False, **kwargs) -> ase.Atoms: 

28 """Create an Atoms instance for a conventional unit cell of a 

29 space group. 

30 

31 Parameters: 

32 

33 symbols : str | sequence of str | sequence of Atom | Atoms 

34 Element symbols of the unique sites. Can either be a string 

35 formula or a sequence of element symbols. E.g. ('Na', 'Cl') 

36 and 'NaCl' are equivalent. Can also be given as a sequence of 

37 Atom objects or an Atoms object. 

38 basis : list of scaled coordinates 

39 Positions of the unique sites corresponding to symbols given 

40 either as scaled positions or through an atoms instance. Not 

41 needed if *symbols* is a sequence of Atom objects or an Atoms 

42 object. 

43 occupancies : list of site occupancies 

44 Occupancies of the unique sites. Defaults to 1.0 and thus no mixed 

45 occupancies are considered if not explicitly asked for. If occupancies 

46 are given, the most dominant species will yield the atomic number. 

47 The occupancies in the atoms.info['occupancy'] dictionary will have 

48 integers keys converted to strings. The conversion is done in order 

49 to avoid unexpected conversions when using the JSON serializer. 

50 spacegroup : int | string | Spacegroup instance 

51 Space group given either as its number in International Tables 

52 or as its Hermann-Mauguin symbol. 

53 setting : 1 | 2 

54 Space group setting. 

55 cell : 3x3 matrix 

56 Unit cell vectors. 

57 cellpar : [a, b, c, alpha, beta, gamma] 

58 Cell parameters with angles in degree. Is not used when `cell` 

59 is given. 

60 ab_normal : vector 

61 Is used to define the orientation of the unit cell relative 

62 to the Cartesian system when `cell` is not given. It is the 

63 normal vector of the plane spanned by a and b. 

64 a_direction : vector 

65 Defines the orientation of the unit cell a vector. a will be 

66 parallel to the projection of `a_direction` onto the a-b plane. 

67 size : 3 positive integers 

68 How many times the conventional unit cell should be repeated 

69 in each direction. 

70 onduplicates : 'keep' | 'replace' | 'warn' | 'error' 

71 Action if `basis` contain symmetry-equivalent positions: 

72 'keep' - ignore additional symmetry-equivalent positions 

73 'replace' - replace 

74 'warn' - like 'keep', but issue an UserWarning 

75 'error' - raises a SpacegroupValueError 

76 symprec : float 

77 Minimum "distance" betweed two sites in scaled coordinates 

78 before they are counted as the same site. 

79 pbc : one or three bools 

80 Periodic boundary conditions flags. Examples: True, 

81 False, 0, 1, (1, 1, 0), (True, False, False). Default 

82 is True. 

83 primitive_cell : bool 

84 Whether to return the primitive instead of the conventional 

85 unit cell. 

86 

87 Keyword arguments: 

88 

89 All additional keyword arguments are passed on to the Atoms 

90 constructor. Currently, probably the most useful additional 

91 keyword arguments are `info`, `constraint` and `calculator`. 

92 

93 Examples: 

94 

95 Two diamond unit cells (space group number 227) 

96 

97 >>> diamond = crystal('C', [(0,0,0)], spacegroup=227, 

98 ... cellpar=[3.57, 3.57, 3.57, 90, 90, 90], size=(2,1,1)) 

99 >>> ase.view(diamond) # doctest: +SKIP 

100 

101 A CoSb3 skutterudite unit cell containing 32 atoms 

102 

103 >>> skutterudite = crystal(('Co', 'Sb'), 

104 ... basis=[(0.25,0.25,0.25), (0.0, 0.335, 0.158)], 

105 ... spacegroup=204, cellpar=[9.04, 9.04, 9.04, 90, 90, 90]) 

106 >>> len(skutterudite) 

107 32 

108 """ 

109 sg = Spacegroup(spacegroup, setting) 

110 if (not isinstance(symbols, str) and 

111 hasattr(symbols, '__getitem__') and 

112 len(symbols) > 0 and 

113 isinstance(symbols[0], ase.Atom)): 

114 symbols = ase.Atoms(symbols) 

115 if isinstance(symbols, ase.Atoms): 

116 basis = symbols 

117 symbols = basis.get_chemical_symbols() 

118 if isinstance(basis, ase.Atoms): 

119 basis_coords = basis.get_scaled_positions() 

120 if cell is None and cellpar is None: 

121 cell = basis.cell 

122 if symbols is None: 

123 symbols = basis.get_chemical_symbols() 

124 else: 

125 basis_coords = np.array(basis, dtype=float, copy=False, ndmin=2) 

126 

127 if occupancies is not None: 

128 occupancies_dict = {} 

129 

130 for index, coord in enumerate(basis_coords): 

131 # Compute all distances and get indices of nearest atoms 

132 dist = spatial.distance.cdist(coord.reshape(1, 3), basis_coords) 

133 indices_dist = np.flatnonzero(dist < symprec) 

134 

135 occ = {symbols[index]: occupancies[index]} 

136 

137 # Check nearest and update occupancy 

138 for index_dist in indices_dist: 

139 if index == index_dist: 

140 continue 

141 else: 

142 occ.update({symbols[index_dist]: occupancies[index_dist]}) 

143 

144 occupancies_dict[str(index)] = occ.copy() 

145 

146 sites, kinds = sg.equivalent_sites(basis_coords, 

147 onduplicates=onduplicates, 

148 symprec=symprec) 

149 

150 # this is needed to handle deuterium masses 

151 masses = None 

152 if 'masses' in kwargs: 

153 masses = kwargs['masses'][kinds] 

154 del kwargs['masses'] 

155 

156 symbols = parse_symbols(symbols) 

157 

158 if occupancies is None: 

159 symbols = [symbols[i] for i in kinds] 

160 else: 

161 # make sure that we put the dominant species there 

162 symbols = [sorted(occupancies_dict[str(i)].items(), key=lambda x: x[1])[-1][0] for i in kinds] 

163 

164 if cell is None: 

165 cell = cellpar_to_cell(cellpar, ab_normal, a_direction) 

166 

167 info: Dict[str, Any] = {} 

168 info['spacegroup'] = sg 

169 if primitive_cell: 

170 info['unit_cell'] = 'primitive' 

171 else: 

172 info['unit_cell'] = 'conventional' 

173 

174 if 'info' in kwargs: 

175 info.update(kwargs['info']) 

176 

177 if occupancies is not None: 

178 info['occupancy'] = occupancies_dict 

179 

180 kwargs['info'] = info 

181 

182 atoms = ase.Atoms(symbols, 

183 scaled_positions=sites, 

184 cell=cell, 

185 pbc=pbc, 

186 masses=masses, 

187 **kwargs) 

188 

189 if isinstance(basis, ase.Atoms): 

190 for name in basis.arrays: 

191 if not atoms.has(name): 

192 array = basis.get_array(name) 

193 atoms.new_array(name, [array[i] for i in kinds], 

194 dtype=array.dtype, shape=array.shape[1:]) 

195 

196 if kinds: 

197 atoms.new_array('spacegroup_kinds', np.asarray(kinds, dtype=int)) 

198 

199 if primitive_cell: 

200 from ase.build import cut 

201 prim_cell = sg.scaled_primitive_cell 

202 

203 # Preserve calculator if present: 

204 calc = atoms.calc 

205 atoms = cut(atoms, a=prim_cell[0], b=prim_cell[1], c=prim_cell[2]) 

206 atoms.calc = calc 

207 

208 if size != (1, 1, 1): 

209 atoms = atoms.repeat(size) 

210 return atoms 

211 

212 

213def parse_symbols(symbols): 

214 """Return `sumbols` as a sequence of element symbols.""" 

215 if isinstance(symbols, str): 

216 symbols = string2symbols(symbols) 

217 return symbols