1"""Helper functions for read_fdf.""" 

2from pathlib import Path 

3from re import compile 

4 

5import numpy as np 

6 

7from ase import Atoms 

8from ase.utils import reader 

9from ase.units import Bohr 

10 

11 

12_label_strip_re = compile(r'[\s._-]') 

13 

14 

15def _labelize(raw_label): 

16 # Labels are case insensitive and -_. should be ignored, lower and strip it 

17 return _label_strip_re.sub('', raw_label).lower() 

18 

19 

20def _is_block(val): 

21 # Tell whether value is a block-value or an ordinary value. 

22 # A block is represented as a list of lists of strings, 

23 # and a ordinary value is represented as a list of strings 

24 if isinstance(val, list) and \ 

25 len(val) > 0 and \ 

26 isinstance(val[0], list): 

27 return True 

28 return False 

29 

30 

31def _get_stripped_lines(fd): 

32 # Remove comments, leading blanks, and empty lines 

33 return [_f for _f in [L.split('#')[0].strip() for L in fd] if _f] 

34 

35 

36@reader 

37def _read_fdf_lines(file): 

38 # Read lines and resolve includes 

39 lbz = _labelize 

40 

41 lines = [] 

42 for L in _get_stripped_lines(file): 

43 w0 = lbz(L.split(None, 1)[0]) 

44 

45 if w0 == '%include': 

46 # Include the contents of fname 

47 fname = L.split(None, 1)[1].strip() 

48 parent_fname = getattr(file, 'name', None) 

49 if isinstance(parent_fname, str): 

50 fname = Path(parent_fname).parent / fname 

51 lines += _read_fdf_lines(fname) 

52 

53 elif '<' in L: 

54 L, fname = L.split('<', 1) 

55 w = L.split() 

56 fname = fname.strip() 

57 

58 if w0 == '%block': 

59 # "%block label < filename" means that the block contents 

60 # should be read from filename 

61 if len(w) != 2: 

62 raise IOError('Bad %%block-statement "%s < %s"' % 

63 (L, fname)) 

64 label = lbz(w[1]) 

65 lines.append('%%block %s' % label) 

66 lines += _get_stripped_lines(open(fname)) 

67 lines.append('%%endblock %s' % label) 

68 else: 

69 # "label < filename.fdf" means that the label 

70 # (_only_ that label) is to be resolved from filename.fdf 

71 label = lbz(w[0]) 

72 fdf = read_fdf(fname) 

73 if label in fdf: 

74 if _is_block(fdf[label]): 

75 lines.append('%%block %s' % label) 

76 lines += [' '.join(x) for x in fdf[label]] 

77 lines.append('%%endblock %s' % label) 

78 else: 

79 lines.append('%s %s' % (label, ' '.join(fdf[label]))) 

80 # else: 

81 # label unresolved! 

82 # One should possibly issue a warning about this! 

83 else: 

84 # Simple include line L 

85 lines.append(L) 

86 return lines 

87 

88 

89def read_fdf(fname): 

90 """Read a siesta style fdf-file. 

91 

92 The data is returned as a dictionary 

93 ( label:value ). 

94 

95 All labels are converted to lower case characters and 

96 are stripped of any '-', '_', or '.'. 

97 

98 Ordinary values are stored as a list of strings (splitted on WS), 

99 and block values are stored as list of lists of strings 

100 (splitted per line, and on WS). 

101 If a label occurres more than once, the first occurrence 

102 takes precedence. 

103 

104 The implementation applies no intelligence, and does not 

105 "understand" the data or the concept of units etc. 

106 Values are never parsed in any way, just stored as 

107 split strings. 

108 

109 The implementation tries to comply with the fdf-format 

110 specification as presented in the siesta 2.0.2 manual. 

111 

112 An fdf-dictionary could e.g. look like this:: 

113 

114 {'atomiccoordinatesandatomicspecies': [ 

115 ['4.9999998', '5.7632392', '5.6095972', '1'], 

116 ['5.0000000', '6.5518100', '4.9929091', '2'], 

117 ['5.0000000', '4.9746683', '4.9929095', '2']], 

118 'atomiccoordinatesformat': ['Ang'], 

119 'chemicalspecieslabel': [['1', '8', 'O'], 

120 ['2', '1', 'H']], 

121 'dmmixingweight': ['0.1'], 

122 'dmnumberpulay': ['5'], 

123 'dmusesavedm': ['True'], 

124 'latticeconstant': ['1.000000', 'Ang'], 

125 'latticevectors': [ 

126 ['10.00000000', '0.00000000', '0.00000000'], 

127 ['0.00000000', '11.52647800', '0.00000000'], 

128 ['0.00000000', '0.00000000', '10.59630900']], 

129 'maxscfiterations': ['120'], 

130 'meshcutoff': ['2721.139566', 'eV'], 

131 'numberofatoms': ['3'], 

132 'numberofspecies': ['2'], 

133 'paobasissize': ['dz'], 

134 'solutionmethod': ['diagon'], 

135 'systemlabel': ['H2O'], 

136 'wavefunckpoints': [['0.0', '0.0', '0.0']], 

137 'writedenchar': ['T'], 

138 'xcauthors': ['PBE'], 

139 'xcfunctional': ['GGA']} 

140 

141 """ 

142 fdf = {} 

143 lbz = _labelize 

144 lines = _read_fdf_lines(fname) 

145 while lines: 

146 w = lines.pop(0).split(None, 1) 

147 if lbz(w[0]) == '%block': 

148 # Block value 

149 if len(w) == 2: 

150 label = lbz(w[1]) 

151 content = [] 

152 while True: 

153 if len(lines) == 0: 

154 raise IOError('Unexpected EOF reached in %s, ' 

155 'un-ended block %s' % (fname, label)) 

156 w = lines.pop(0).split() 

157 if lbz(w[0]) == '%endblock': 

158 break 

159 content.append(w) 

160 

161 if label not in fdf: 

162 # Only first appearance of label is to be used 

163 fdf[label] = content 

164 else: 

165 raise IOError('%%block statement without label') 

166 else: 

167 # Ordinary value 

168 label = lbz(w[0]) 

169 if len(w) == 1: 

170 # Siesta interpret blanks as True for logical variables 

171 fdf[label] = [] 

172 else: 

173 fdf[label] = w[1].split() 

174 return fdf 

175 

176 

177def read_struct_out(fd): 

178 """Read a siesta struct file""" 

179 

180 cell = [] 

181 for i in range(3): 

182 line = next(fd) 

183 v = np.array(line.split(), float) 

184 cell.append(v) 

185 

186 natoms = int(next(fd)) 

187 

188 numbers = np.empty(natoms, int) 

189 scaled_positions = np.empty((natoms, 3)) 

190 for i, line in enumerate(fd): 

191 tokens = line.split() 

192 numbers[i] = int(tokens[1]) 

193 scaled_positions[i] = np.array(tokens[2:5], float) 

194 

195 return Atoms(numbers, 

196 cell=cell, 

197 pbc=True, 

198 scaled_positions=scaled_positions) 

199 

200 

201def read_siesta_xv(fd): 

202 vectors = [] 

203 for i in range(3): 

204 data = next(fd).split() 

205 vectors.append([float(data[j]) * Bohr for j in range(3)]) 

206 

207 # Read number of atoms (line 4) 

208 natoms = int(next(fd).split()[0]) 

209 

210 # Read remaining lines 

211 speciesnumber, atomnumbers, xyz, V = [], [], [], [] 

212 for line in fd: 

213 if len(line) > 5: # Ignore blank lines 

214 data = line.split() 

215 speciesnumber.append(int(data[0])) 

216 atomnumbers.append(int(data[1])) 

217 xyz.append([float(data[2 + j]) * Bohr for j in range(3)]) 

218 V.append([float(data[5 + j]) * Bohr for j in range(3)]) 

219 

220 vectors = np.array(vectors) 

221 atomnumbers = np.array(atomnumbers) 

222 xyz = np.array(xyz) 

223 atoms = Atoms(numbers=atomnumbers, positions=xyz, cell=vectors, 

224 pbc=True) 

225 assert natoms == len(atoms) 

226 return atoms