Coverage for /builds/ase/ase/ase/io/siesta.py: 72.90%

1# fmt: off 

2 

3"""Helper functions for read_fdf.""" 

4from pathlib import Path 

5from re import compile 

6 

7import numpy as np 

8 

9from ase import Atoms 

10from ase.units import Bohr 

11from ase.utils import reader 

12 

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

14 

15 

16def _labelize(raw_label): 

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

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

19 

20 

21def _is_block(val): 

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

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

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

25 if isinstance(val, list) and \ 

26 len(val) > 0 and \ 

27 isinstance(val[0], list): 

28 return True 

29 return False 

30 

31 

32def _get_stripped_lines(fd): 

33 # Remove comments, leading blanks, and empty lines 

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

35 

36 

37@reader 

38def _read_fdf_lines(file): 

39 # Read lines and resolve includes 

40 lbz = _labelize 

41 

42 lines = [] 

43 for L in _get_stripped_lines(file): 

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

45 

46 if w0 == '%include': 

47 # Include the contents of fname 

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

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

50 if isinstance(parent_fname, str): 

51 fname = Path(parent_fname).parent / fname 

52 lines += _read_fdf_lines(fname) 

53 

54 elif '<' in L: 

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

56 w = L.split() 

57 fname = fname.strip() 

58 

59 if w0 == '%block': 

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

61 # should be read from filename 

62 if len(w) != 2: 

63 raise OSError('Bad %%block-statement "%s < %s"' % 

64 (L, fname)) 

65 label = lbz(w[1]) 

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

67 lines += _get_stripped_lines(open(fname)) 

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

69 else: 

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

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

72 label = lbz(w[0]) 

73 fdf = read_fdf(fname) 

74 if label in fdf: 

75 if _is_block(fdf[label]): 

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

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

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

79 else: 

80 lines.append('{} {}'.format( 

81 label, ' '.join(fdf[label]))) 

82 # else: 

83 # label unresolved! 

84 # One should possibly issue a warning about this! 

85 else: 

86 # Simple include line L 

87 lines.append(L) 

88 return lines 

89 

90 

91def read_fdf(fname): 

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

93 

94 The data is returned as a dictionary 

95 ( label:value ). 

96 

97 All labels are converted to lower case characters and 

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

99 

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

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

102 (splitted per line, and on WS). 

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

104 takes precedence. 

105 

106 The implementation applies no intelligence, and does not 

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

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

109 split strings. 

110 

111 The implementation tries to comply with the fdf-format 

112 specification as presented in the siesta 2.0.2 manual. 

113 

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

115 

116 {'atomiccoordinatesandatomicspecies': [ 

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

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

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

120 'atomiccoordinatesformat': ['Ang'], 

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

122 ['2', '1', 'H']], 

123 'dmmixingweight': ['0.1'], 

124 'dmnumberpulay': ['5'], 

125 'dmusesavedm': ['True'], 

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

127 'latticevectors': [ 

128 ['10.00000000', '0.00000000', '0.00000000'], 

129 ['0.00000000', '11.52647800', '0.00000000'], 

130 ['0.00000000', '0.00000000', '10.59630900']], 

131 'maxscfiterations': ['120'], 

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

133 'numberofatoms': ['3'], 

134 'numberofspecies': ['2'], 

135 'paobasissize': ['dz'], 

136 'solutionmethod': ['diagon'], 

137 'systemlabel': ['H2O'], 

138 'wavefunckpoints': [['0.0', '0.0', '0.0']], 

139 'writedenchar': ['T'], 

140 'xcauthors': ['PBE'], 

141 'xcfunctional': ['GGA']} 

142 

143 """ 

144 fdf = {} 

145 lbz = _labelize 

146 lines = _read_fdf_lines(fname) 

147 while lines: 

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

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

150 # Block value 

151 if len(w) == 2: 

152 label = lbz(w[1]) 

153 content = [] 

154 while True: 

155 if len(lines) == 0: 

156 raise OSError('Unexpected EOF reached in %s, ' 

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

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

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

160 break 

161 content.append(w) 

162 

163 if label not in fdf: 

164 # Only first appearance of label is to be used 

165 fdf[label] = content 

166 else: 

167 raise OSError('%%block statement without label') 

168 else: 

169 # Ordinary value 

170 label = lbz(w[0]) 

171 if len(w) == 1: 

172 # Siesta interpret blanks as True for logical variables 

173 fdf[label] = [] 

174 else: 

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

176 return fdf 

177 

178 

179def read_struct_out(fd): 

180 """Read a siesta struct file""" 

181 

182 cell = [] 

183 for _ in range(3): 

184 line = next(fd) 

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

186 cell.append(v) 

187 

188 natoms = int(next(fd)) 

189 

190 numbers = np.empty(natoms, int) 

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

192 for i, line in enumerate(fd): 

193 tokens = line.split() 

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

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

196 

197 return Atoms(numbers, 

198 cell=cell, 

199 pbc=True, 

200 scaled_positions=scaled_positions) 

201 

202 

203def read_siesta_xv(fd): 

204 vectors = [] 

205 for _ in range(3): 

206 data = next(fd).split() 

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

208 

209 # Read number of atoms (line 4) 

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

211 

212 # Read remaining lines 

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

214 for line in fd: 

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

216 data = line.split() 

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

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

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

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

221 

222 vectors = np.array(vectors) 

223 atomnumbers = np.array(atomnumbers) 

224 xyz = np.array(xyz) 

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

226 pbc=True) 

227 assert natoms == len(atoms) 

228 return atoms