Coverage for /builds/ase/ase/ase/spacegroup/xtal.py: 93.33%

1# fmt: off 

2 

3# Copyright (C) 2010, Jesper Friis 

4# (see accompanying license files for details). 

5 

6""" 

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

8knowledge of the space group. 

9 

10""" 

11 

12from typing import Any, Dict 

13 

14import numpy as np 

15from scipy import spatial 

16 

17import ase 

18from ase.geometry import cellpar_to_cell 

19from ase.spacegroup import Spacegroup 

20from ase.symbols import string2symbols 

21 

22__all__ = ['crystal'] 

23 

24 

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

26 cell=None, cellpar=None, 

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

28 onduplicates='warn', symprec=0.001, 

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

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

31 space group. 

32 

33 Parameters: 

34 

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

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

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

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

39 Atom objects or an Atoms object. 

40 basis : list of scaled coordinates 

41 Positions of the unique sites corresponding to symbols given 

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

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

44 object. 

45 occupancies : list of site occupancies 

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

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

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

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

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

51 to avoid unexpected conversions when using the JSON serializer. 

52 spacegroup : int | string | Spacegroup instance 

53 Space group given either as its number in International Tables 

54 or as its Hermann-Mauguin symbol. 

55 setting : 1 | 2 

56 Space group setting. 

57 cell : 3x3 matrix 

58 Unit cell vectors. 

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

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

61 is given. 

62 ab_normal : vector 

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

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

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

66 a_direction : vector 

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

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

69 size : 3 positive integers 

70 How many times the conventional unit cell should be repeated 

71 in each direction. 

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

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

74 'keep' - ignore additional symmetry-equivalent positions 

75 'replace' - replace 

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

77 'error' - raises a SpacegroupValueError 

78 symprec : float 

79 Minimum "distance" betweed two sites in scaled coordinates 

80 before they are counted as the same site. 

81 pbc : one or three bools 

82 Periodic boundary conditions flags. Examples: True, 

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

84 is True. 

85 primitive_cell : bool 

86 Whether to return the primitive instead of the conventional 

87 unit cell. 

88 

89 Keyword arguments: 

90 

91 All additional keyword arguments are passed on to the Atoms 

92 constructor. Currently, probably the most useful additional 

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

94 

95 Examples: 

96 

97 Two diamond unit cells (space group number 227) 

98 

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

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

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

102 

103 A CoSb3 skutterudite unit cell containing 32 atoms 

104 

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

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

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

108 >>> len(skutterudite) 

109 32 

110 """ 

111 sg = Spacegroup(spacegroup, setting) 

112 if ( 

113 not isinstance(symbols, str) and 

114 hasattr(symbols, '__getitem__') and 

115 len(symbols) > 0 and 

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

117 symbols = ase.Atoms(symbols) 

118 if isinstance(symbols, ase.Atoms): 

119 basis = symbols 

120 symbols = basis.get_chemical_symbols() 

121 if isinstance(basis, ase.Atoms): 

122 basis_coords = basis.get_scaled_positions() 

123 if cell is None and cellpar is None: 

124 cell = basis.cell 

125 if symbols is None: 

126 symbols = basis.get_chemical_symbols() 

127 else: 

128 basis_coords = np.array(basis, dtype=float, ndmin=2) 

129 

130 if occupancies is not None: 

131 occupancies_dict = {} 

132 

133 for index, coord in enumerate(basis_coords): 

134 # Compute all distances and get indices of nearest atoms 

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

136 indices_dist = np.flatnonzero(dist < symprec) 

137 

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

139 

140 # Check nearest and update occupancy 

141 for index_dist in indices_dist: 

142 if index == index_dist: 

143 continue 

144 else: 

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

146 

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

148 

149 sites, kinds = sg.equivalent_sites(basis_coords, 

150 onduplicates=onduplicates, 

151 symprec=symprec) 

152 

153 # this is needed to handle deuterium masses 

154 masses = None 

155 if 'masses' in kwargs: 

156 masses = kwargs['masses'][kinds] 

157 del kwargs['masses'] 

158 

159 symbols = parse_symbols(symbols) 

160 

161 if occupancies is None: 

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

163 else: 

164 # make sure that we put the dominant species there 

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

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

167 

168 if cell is None: 

169 cell = cellpar_to_cell(cellpar, ab_normal, a_direction) 

170 

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

172 info['spacegroup'] = sg 

173 if primitive_cell: 

174 info['unit_cell'] = 'primitive' 

175 else: 

176 info['unit_cell'] = 'conventional' 

177 

178 if 'info' in kwargs: 

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

180 

181 if occupancies is not None: 

182 info['occupancy'] = occupancies_dict 

183 

184 kwargs['info'] = info 

185 

186 atoms = ase.Atoms(symbols, 

187 scaled_positions=sites, 

188 cell=cell, 

189 pbc=pbc, 

190 masses=masses, 

191 **kwargs) 

192 

193 if isinstance(basis, ase.Atoms): 

194 for name in basis.arrays: 

195 if not atoms.has(name): 

196 array = basis.get_array(name) 

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

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

199 

200 if kinds: 

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

202 

203 if primitive_cell: 

204 from ase.build import cut 

205 prim_cell = sg.scaled_primitive_cell 

206 

207 # Preserve calculator if present: 

208 calc = atoms.calc 

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

210 atoms.calc = calc 

211 

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

213 atoms = atoms.repeat(size) 

214 return atoms 

215 

216 

217def parse_symbols(symbols): 

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

219 if isinstance(symbols, str): 

220 symbols = string2symbols(symbols) 

221 return symbols