Coverage for /builds/kinetik161/ase/ase/spacegroup/xtal.py: 93.42%

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 Any, Dict 

11 

12import numpy as np 

13from scipy import spatial 

14 

15import ase 

16from ase.geometry import cellpar_to_cell 

17from ase.spacegroup import Spacegroup 

18from ase.symbols import string2symbols 

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 ( 

111 not isinstance(symbols, str) and 

112 hasattr(symbols, '__getitem__') and 

113 len(symbols) > 0 and 

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

115 symbols = ase.Atoms(symbols) 

116 if isinstance(symbols, ase.Atoms): 

117 basis = symbols 

118 symbols = basis.get_chemical_symbols() 

119 if isinstance(basis, ase.Atoms): 

120 basis_coords = basis.get_scaled_positions() 

121 if cell is None and cellpar is None: 

122 cell = basis.cell 

123 if symbols is None: 

124 symbols = basis.get_chemical_symbols() 

125 else: 

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

127 

128 if occupancies is not None: 

129 occupancies_dict = {} 

130 

131 for index, coord in enumerate(basis_coords): 

132 # Compute all distances and get indices of nearest atoms 

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

134 indices_dist = np.flatnonzero(dist < symprec) 

135 

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

137 

138 # Check nearest and update occupancy 

139 for index_dist in indices_dist: 

140 if index == index_dist: 

141 continue 

142 else: 

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

144 

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

146 

147 sites, kinds = sg.equivalent_sites(basis_coords, 

148 onduplicates=onduplicates, 

149 symprec=symprec) 

150 

151 # this is needed to handle deuterium masses 

152 masses = None 

153 if 'masses' in kwargs: 

154 masses = kwargs['masses'][kinds] 

155 del kwargs['masses'] 

156 

157 symbols = parse_symbols(symbols) 

158 

159 if occupancies is None: 

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

161 else: 

162 # make sure that we put the dominant species there 

163 symbols = [sorted(occupancies_dict[str(i)].items(), 

164 key=lambda x: x[1])[-1][0] for i in kinds] 

165 

166 if cell is None: 

167 cell = cellpar_to_cell(cellpar, ab_normal, a_direction) 

168 

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

170 info['spacegroup'] = sg 

171 if primitive_cell: 

172 info['unit_cell'] = 'primitive' 

173 else: 

174 info['unit_cell'] = 'conventional' 

175 

176 if 'info' in kwargs: 

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

178 

179 if occupancies is not None: 

180 info['occupancy'] = occupancies_dict 

181 

182 kwargs['info'] = info 

183 

184 atoms = ase.Atoms(symbols, 

185 scaled_positions=sites, 

186 cell=cell, 

187 pbc=pbc, 

188 masses=masses, 

189 **kwargs) 

190 

191 if isinstance(basis, ase.Atoms): 

192 for name in basis.arrays: 

193 if not atoms.has(name): 

194 array = basis.get_array(name) 

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

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

197 

198 if kinds: 

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

200 

201 if primitive_cell: 

202 from ase.build import cut 

203 prim_cell = sg.scaled_primitive_cell 

204 

205 # Preserve calculator if present: 

206 calc = atoms.calc 

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

208 atoms.calc = calc 

209 

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

211 atoms = atoms.repeat(size) 

212 return atoms 

213 

214 

215def parse_symbols(symbols): 

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

217 if isinstance(symbols, str): 

218 symbols = string2symbols(symbols) 

219 return symbols