Coverage for /builds/debichem-team/python-ase/ase/calculators/openmx/reader.py: 63.66%

1# flake8: noqa 

2""" 

3The ASE Calculator for OpenMX <http://www.openmx-square.org>: Python interface 

4to the software package for nano-scale material simulations based on density 

5functional theories. 

6 Copyright (C) 2018 JaeHwan Shim and JaeJun Yu 

7 

8 This program is free software: you can redistribute it and/or modify 

9 it under the terms of the GNU Lesser General Public License as published by 

10 the Free Software Foundation, either version 2.1 of the License, or 

11 (at your option) any later version. 

12 

13 This program is distributed in the hope that it will be useful, 

14 but WITHOUT ANY WARRANTY; without even the implied warranty of 

15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

16 GNU Lesser General Public License for more details. 

17 

18 You should have received a copy of the GNU Lesser General Public License 

19 along with ASE. If not, see <http://www.gnu.org/licenses/>. 

20""" 

21# from ase.calculators import SinglePointDFTCalculator 

22import os 

23import struct 

24 

25import numpy as np 

26 

27from ase.io import ParseError 

28from ase.units import Bohr, Debye, Ha 

29 

30 

31def read_openmx(filename=None, debug=False): 

32 from ase import Atoms 

33 from ase.calculators.openmx import OpenMX 

34 """ 

35 Read results from typical OpenMX output files and returns the atom object 

36 In default mode, it reads every implementd properties we could get from 

37 the files. Unlike previous version, we read the information based on file. 

38 previous results will be eraised unless previous results are written in the 

39 next calculation results. 

40 

41 Read the 'LABEL.log' file seems redundant. Because all the 

42 information should already be written in '.out' file. However, in the 

43 version 3.8.3, stress tensor are not written in the '.out' file. It only 

44 contained in the '.log' file. So... I implented reading '.log' file method 

45 """ 

46 log_data = read_file(get_file_name('.log', filename), debug=debug) 

47 restart_data = read_file(get_file_name('.dat#', filename), debug=debug) 

48 dat_data = read_file(get_file_name('.dat', filename), debug=debug) 

49 out_data = read_file(get_file_name('.out', filename), debug=debug) 

50 scfout_data = read_scfout_file(get_file_name('.scfout', filename)) 

51 band_data = read_band_file(get_file_name('.Band', filename)) 

52 # dos_data = read_dos_file(get_file_name('.Dos.val', filename)) 

53 """ 

54 First, we get every data we could get from the all results files. And then, 

55 reform the data to fit to data structure of Atom object. While doing this, 

56 Fix the unit to ASE format. 

57 """ 

58 parameters = get_parameters(out_data=out_data, log_data=log_data, 

59 restart_data=restart_data, dat_data=dat_data, 

60 scfout_data=scfout_data, band_data=band_data) 

61 atomic_formula = get_atomic_formula(out_data=out_data, log_data=log_data, 

62 restart_data=restart_data, 

63 scfout_data=scfout_data, 

64 dat_data=dat_data) 

65 results = get_results(out_data=out_data, log_data=log_data, 

66 restart_data=restart_data, scfout_data=scfout_data, 

67 dat_data=dat_data, band_data=band_data) 

68 

69 atoms = Atoms(**atomic_formula) 

70 

71 # XXX Creating OpenMX out of parameters directly is a security problem 

72 # since the data comes from files, and the files may contain 

73 # malicious content that would be interpreted as 'command' 

74 from ase.calculators.calculator import all_properties 

75 from ase.calculators.singlepoint import SinglePointDFTCalculator 

76 common_properties = set(all_properties) & set(results) 

77 results = {key: results[key] for key in common_properties} 

78 

79 atoms.calc = SinglePointDFTCalculator( 

80 atoms, **results) 

81 

82 # atoms.calc = OpenMX(**parameters) 

83 # atoms.calc.results = results 

84 return atoms 

85 

86 

87def read_file(filename, debug=False): 

88 """ 

89 Read the 'LABEL.out' file. Using 'parameters.py', we read every 'allowed_ 

90 dat' dictionory. while reading a file, if one find the key matcheds That 

91 'patters', which indicates the property we want is written, it will returns 

92 the pair value of that key. For example, 

93 example will be written later 

94 """ 

95 from ase.calculators.openmx import parameters as param 

96 if not os.path.isfile(filename): 

97 return {} 

98 param_keys = ['integer_keys', 'float_keys', 'string_keys', 'bool_keys', 

99 'list_int_keys', 'list_float_keys', 'list_bool_keys', 

100 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys'] 

101 patterns = { 

102 'Stress tensor': ('stress', read_stress_tensor), 

103 'Dipole moment': ('dipole', read_dipole), 

104 'Fractional coordinates of': ('scaled_positions', read_scaled_positions), 

105 'Utot.': ('energy', read_energy), 

106 'energies in': ('energies', read_energies), 

107 'Chemical Potential': ('chemical_potential', read_chemical_potential), 

108 '<coordinates.forces': ('forces', read_forces), 

109 'Eigenvalues (Hartree)': ('eigenvalues', read_eigenvalues)} 

110 special_patterns = { 

111 'Total spin moment': (('magmoms', 'total_magmom'), 

112 read_magmoms_and_total_magmom), 

113 } 

114 out_data = {} 

115 line = '\n' 

116 if (debug): 

117 print(f'Read results from {filename}') 

118 with open(filename) as fd: 

119 ''' 

120 Read output file line by line. When the `line` matches the pattern 

121 of certain keywords in `param.[dtype]_keys`, for example, 

122 

123 if line in param.string_keys: 

124 out_data[key] = read_string(line) 

125 

126 parse that line and store it to `out_data` in specified data type. 

127 To cover all `dtype` parameters, for loop was used, 

128 

129 for [dtype] in parameters_keys: 

130 if line in param.[dtype]_keys: 

131 out_data[key] = read_[dtype](line) 

132 

133 After found matched pattern, escape the for loop using `continue`. 

134 ''' 

135 while line != '': 

136 pattern_matched = False 

137 line = fd.readline() 

138 try: 

139 _line = line.split()[0] 

140 except IndexError: 

141 continue 

142 for dtype_key in param_keys: 

143 dtype = dtype_key.rsplit('_', 1)[0] 

144 read_dtype = globals()['read_' + dtype] 

145 for key in param.__dict__[dtype_key]: 

146 if key in _line: 

147 out_data[get_standard_key(key)] = read_dtype(line) 

148 pattern_matched = True 

149 continue 

150 if pattern_matched: 

151 continue 

152 

153 for key in param.matrix_keys: 

154 if '<' + key in line: 

155 out_data[get_standard_key(key)] = read_matrix(line, key, fd) 

156 pattern_matched = True 

157 continue 

158 if pattern_matched: 

159 continue 

160 for key in patterns: 

161 if key in line: 

162 out_data[patterns[key][0]] = patterns[key][1]( 

163 line, fd, debug=debug) 

164 pattern_matched = True 

165 continue 

166 if pattern_matched: 

167 continue 

168 for key in special_patterns: 

169 if key in line: 

170 a, b = special_patterns[key][1](line, fd) 

171 out_data[special_patterns[key][0][0]] = a 

172 out_data[special_patterns[key][0][1]] = b 

173 pattern_matched = True 

174 continue 

175 if pattern_matched: 

176 continue 

177 return out_data 

178 

179 

180def read_scfout_file(filename=None): 

181 """ 

182 Read the Developer output '.scfout' files. It Behaves like read_scfout.c, 

183 OpenMX module, but written in python. Note that some array are begin with 

184 1, not 0 

185 

186 atomnum: the number of total atoms 

187 Catomnum: the number of atoms in the central region 

188 Latomnum: the number of atoms in the left lead 

189 Ratomnum: the number of atoms in the left lead 

190 SpinP_switch: 

191 0: non-spin polarized 

192 1: spin polarized 

193 TCpyCell: the total number of periodic cells 

194 Solver: method for solving eigenvalue problem 

195 ChemP: chemical potential 

196 Valence_Electrons: total number of valence electrons 

197 Total_SpinS: total value of Spin (2*Total_SpinS = muB) 

198 E_Temp: electronic temperature 

199 Total_NumOrbs: the number of atomic orbitals in each atom 

200 size: Total_NumOrbs[atomnum+1] 

201 FNAN: the number of first neighboring atoms of each atom 

202 size: FNAN[atomnum+1] 

203 natn: global index of neighboring atoms of an atom ct_AN 

204 size: natn[atomnum+1][FNAN[ct_AN]+1] 

205 ncn: global index for cell of neighboring atoms of an atom ct_AN 

206 size: ncn[atomnum+1][FNAN[ct_AN]+1] 

207 atv: x,y,and z-components of translation vector of periodically copied cell 

208 size: atv[TCpyCell+1][4]: 

209 atv_ijk: i,j,and j number of periodically copied cells 

210 size: atv_ijk[TCpyCell+1][4]: 

211 tv[4][4]: unit cell vectors in Bohr 

212 rtv[4][4]: reciprocal unit cell vectors in Bohr^{-1} 

213 note: 

214 tv_i dot rtv_j = 2PI * Kronecker's delta_{ij} 

215 Gxyz[atomnum+1][60]: atomic coordinates in Bohr 

216 Hks: Kohn-Sham matrix elements of basis orbitals 

217 size: Hks[SpinP_switch+1] 

218 [atomnum+1] 

219 [FNAN[ct_AN]+1] 

220 [Total_NumOrbs[ct_AN]] 

221 [Total_NumOrbs[h_AN]] 

222 iHks: 

223 imaginary Kohn-Sham matrix elements of basis orbitals 

224 for alpha-alpha, beta-beta, and alpha-beta spin matrices 

225 of which contributions come from spin-orbit coupling 

226 and Hubbard U effective potential. 

227 size: iHks[3] 

228 [atomnum+1] 

229 [FNAN[ct_AN]+1] 

230 [Total_NumOrbs[ct_AN]] 

231 [Total_NumOrbs[h_AN]] 

232 OLP: overlap matrix 

233 size: OLP[atomnum+1] 

234 [FNAN[ct_AN]+1] 

235 [Total_NumOrbs[ct_AN]] 

236 [Total_NumOrbs[h_AN]] 

237 OLPpox: overlap matrix with position operator x 

238 size: OLPpox[atomnum+1] 

239 [FNAN[ct_AN]+1] 

240 [Total_NumOrbs[ct_AN]] 

241 [Total_NumOrbs[h_AN]] 

242 OLPpoy: overlap matrix with position operator y 

243 size: OLPpoy[atomnum+1] 

244 [FNAN[ct_AN]+1] 

245 [Total_NumOrbs[ct_AN]] 

246 [Total_NumOrbs[h_AN]] 

247 OLPpoz: overlap matrix with position operator z 

248 size: OLPpoz[atomnum+1] 

249 [FNAN[ct_AN]+1] 

250 [Total_NumOrbs[ct_AN]] 

251 [Total_NumOrbs[h_AN]] 

252 DM: overlap matrix 

253 size: DM[SpinP_switch+1] 

254 [atomnum+1] 

255 [FNAN[ct_AN]+1] 

256 [Total_NumOrbs[ct_AN]] 

257 [Total_NumOrbs[h_AN]] 

258 dipole_moment_core[4]: 

259 dipole_moment_background[4]: 

260 """ 

261 from numpy import cumsum as cum 

262 from numpy import insert as ins 

263 from numpy import split as spl 

264 from numpy import sum, zeros 

265 if not os.path.isfile(filename): 

266 return {} 

267 

268 def easyReader(byte, data_type, shape): 

269 data_size = {'d': 8, 'i': 4} 

270 data_struct = {'d': float, 'i': int} 

271 dt = data_type 

272 ds = data_size[data_type] 

273 unpack = struct.unpack 

274 if len(byte) == ds: 

275 if dt == 'i': 

276 return data_struct[dt].from_bytes(byte, byteorder='little') 

277 elif dt == 'd': 

278 return np.array(unpack(dt * (len(byte) // ds), byte))[0] 

279 elif shape is not None: 

280 return np.array(unpack(dt * (len(byte) // ds), byte)).reshape(shape) 

281 else: 

282 return np.array(unpack(dt * (len(byte) // ds), byte)) 

283 

284 def inte(byte, shape=None): 

285 return easyReader(byte, 'i', shape) 

286 

287 def floa(byte, shape=None): 

288 return easyReader(byte, 'd', shape) 

289 

290 def readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd): 

291 myOLP = [] 

292 myOLP.append([]) 

293 for ct_AN in range(1, atomnum + 1): 

294 myOLP.append([]) 

295 TNO1 = Total_NumOrbs[ct_AN] 

296 for h_AN in range(FNAN[ct_AN] + 1): 

297 myOLP[ct_AN].append([]) 

298 Gh_AN = natn[ct_AN][h_AN] 

299 TNO2 = Total_NumOrbs[Gh_AN] 

300 for i in range(TNO1): 

301 myOLP[ct_AN][h_AN].append(floa(fd.read(8 * TNO2))) 

302 return myOLP 

303 

304 def readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd): 

305 Hks = [] 

306 for spin in range(SpinP_switch + 1): 

307 Hks.append([]) 

308 Hks[spin].append([np.zeros(FNAN[0] + 1)]) 

309 for ct_AN in range(1, atomnum + 1): 

310 Hks[spin].append([]) 

311 TNO1 = Total_NumOrbs[ct_AN] 

312 for h_AN in range(FNAN[ct_AN] + 1): 

313 Hks[spin][ct_AN].append([]) 

314 Gh_AN = natn[ct_AN][h_AN] 

315 TNO2 = Total_NumOrbs[Gh_AN] 

316 for i in range(TNO1): 

317 Hks[spin][ct_AN][h_AN].append(floa(fd.read(8 * TNO2))) 

318 return Hks 

319 

320 with open(filename, mode='rb') as fd: 

321 atomnum, SpinP_switch = inte(fd.read(8)) 

322 Catomnum, Latomnum, Ratomnum, TCpyCell = inte(fd.read(16)) 

323 atv = floa(fd.read(8 * 4 * (TCpyCell + 1)), shape=(TCpyCell + 1, 4)) 

324 atv_ijk = inte(fd.read(4 * 4 * (TCpyCell + 1)), shape=(TCpyCell + 1, 4)) 

325 Total_NumOrbs = np.insert(inte(fd.read(4 * (atomnum))), 0, 1, axis=0) 

326 FNAN = np.insert(inte(fd.read(4 * (atomnum))), 0, 0, axis=0) 

327 natn = ins(spl(inte(fd.read(4 * sum(FNAN[1:] + 1))), cum(FNAN[1:] + 1)), 

328 0, zeros(FNAN[0] + 1), axis=0)[:-1] 

329 ncn = ins(spl(inte(fd.read(4 * np.sum(FNAN[1:] + 1))), cum(FNAN[1:] + 1)), 

330 0, np.zeros(FNAN[0] + 1), axis=0)[:-1] 

331 tv = ins(floa(fd.read(8 * 3 * 4), shape=(3, 4)), 

332 0, [0, 0, 0, 0], axis=0) 

333 rtv = ins(floa(fd.read(8 * 3 * 4), shape=(3, 4)), 

334 0, [0, 0, 0, 0], axis=0) 

335 Gxyz = ins(floa(fd.read(8 * (atomnum) * 4), shape=(atomnum, 4)), 0, 

336 [0., 0., 0., 0.], axis=0) 

337 Hks = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd) 

338 iHks = [] 

339 if SpinP_switch == 3: 

340 iHks = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd) 

341 OLP = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 

342 OLPpox = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 

343 OLPpoy = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 

344 OLPpoz = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 

345 DM = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd) 

346 Solver = inte(fd.read(4)) 

347 ChemP, E_Temp = floa(fd.read(8 * 2)) 

348 dipole_moment_core = floa(fd.read(8 * 3)) 

349 dipole_moment_background = floa(fd.read(8 * 3)) 

350 Valence_Electrons, Total_SpinS = floa(fd.read(8 * 2)) 

351 

352 scf_out = {'atomnum': atomnum, 'SpinP_switch': SpinP_switch, 

353 'Catomnum': Catomnum, 'Latomnum': Latomnum, 'Hks': Hks, 

354 'Ratomnum': Ratomnum, 'TCpyCell': TCpyCell, 'atv': atv, 

355 'Total_NumOrbs': Total_NumOrbs, 'FNAN': FNAN, 'natn': natn, 

356 'ncn': ncn, 'tv': tv, 'rtv': rtv, 'Gxyz': Gxyz, 'OLP': OLP, 

357 'OLPpox': OLPpox, 'OLPpoy': OLPpoy, 'OLPpoz': OLPpoz, 

358 'Solver': Solver, 'ChemP': ChemP, 'E_Temp': E_Temp, 

359 'dipole_moment_core': dipole_moment_core, 'iHks': iHks, 

360 'dipole_moment_background': dipole_moment_background, 

361 'Valence_Electrons': Valence_Electrons, 'atv_ijk': atv_ijk, 

362 'Total_SpinS': Total_SpinS, 'DM': DM 

363 } 

364 return scf_out 

365 

366 

367def read_band_file(filename=None): 

368 band_data = {} 

369 if not os.path.isfile(filename): 

370 return {} 

371 band_kpath = [] 

372 eigen_bands = [] 

373 with open(filename) as fd: 

374 line = fd.readline().split() 

375 nkpts = 0 

376 nband = int(line[0]) 

377 nspin = int(line[1]) + 1 

378 band_data['nband'] = nband 

379 band_data['nspin'] = nspin 

380 line = fd.readline().split() 

381 band_data['band_kpath_unitcell'] = [line[:3], line[3:6], line[6:9]] 

382 line = fd.readline().split() 

383 band_data['band_nkpath'] = int(line[0]) 

384 for _ in range(band_data['band_nkpath']): 

385 line = fd.readline().split() 

386 band_kpath.append(line) 

387 nkpts += int(line[0]) 

388 band_data['nkpts'] = nkpts 

389 band_data['band_kpath'] = band_kpath 

390 kpts = np.zeros((nkpts, 3)) 

391 eigen_bands = np.zeros((nspin, nkpts, nband)) 

392 for i in range(nspin): 

393 for j in range(nkpts): 

394 line = fd.readline() 

395 kpts[j] = np.array(line.split(), dtype=float)[1:] 

396 line = fd.readline() 

397 eigen_bands[i, j] = np.array(line.split(), dtype=float)[:] 

398 band_data['eigenvalues'] = eigen_bands 

399 band_data['band_kpts'] = kpts 

400 return band_data 

401 

402 

403def read_electron_valency(filename='H_CA13'): 

404 array = [] 

405 with open(filename) as fd: 

406 array = fd.readlines() 

407 fd.close() 

408 required_line = '' 

409 for line in array: 

410 if 'valence.electron' in line: 

411 required_line = line 

412 return rn(required_line) 

413 

414 

415def rn(line='\n', n=1): 

416 """ 

417 Read n'th to last value. 

418 For example: 

419 ... 

420 scf.XcType LDA 

421 scf.Kgrid 4 4 4 

422 ... 

423 In Python, 

424 >>> str(rn(line, 1)) 

425 LDA 

426 >>> line = fd.readline() 

427 >>> int(rn(line, 3)) 

428 4 

429 """ 

430 return line.split()[-n] 

431 

432 

433def read_tuple_integer(line): 

434 return tuple(int(x) for x in line.split()[-3:]) 

435 

436 

437def read_tuple_float(line): 

438 return tuple(float(x) for x in line.split()[-3:]) 

439 

440 

441def read_integer(line): 

442 return int(rn(line)) 

443 

444 

445def read_float(line): 

446 return float(rn(line)) 

447 

448 

449def read_string(line): 

450 return str(rn(line)) 

451 

452 

453def read_bool(line): 

454 bool = str(rn(line)).lower() 

455 if bool == 'on': 

456 return True 

457 elif bool == 'off': 

458 return False 

459 else: 

460 return None 

461 

462 

463def read_list_int(line): 

464 return [int(x) for x in line.split()[1:]] 

465 

466 

467def read_list_float(line): 

468 return [float(x) for x in line.split()[1:]] 

469 

470 

471def read_list_bool(line): 

472 return [read_bool(x) for x in line.split()[1:]] 

473 

474 

475def read_matrix(line, key, fd): 

476 matrix = [] 

477 line = fd.readline() 

478 while key not in line: 

479 matrix.append(line.split()) 

480 line = fd.readline() 

481 return matrix 

482 

483 

484def read_stress_tensor(line, fd, debug=None): 

485 fd.readline() # passing empty line 

486 fd.readline() 

487 line = fd.readline() 

488 xx, xy, xz = read_tuple_float(line) 

489 line = fd.readline() 

490 yx, yy, yz = read_tuple_float(line) 

491 line = fd.readline() 

492 zx, zy, zz = read_tuple_float(line) 

493 stress = [xx, yy, zz, (zy + yz) / 2, (zx + xz) / 2, (yx + xy) / 2] 

494 return stress 

495 

496 

497def read_magmoms_and_total_magmom(line, fd, debug=None): 

498 total_magmom = read_float(line) 

499 fd.readline() # Skip empty lines 

500 fd.readline() 

501 line = fd.readline() 

502 magmoms = [] 

503 while not (line == '' or line.isspace()): 

504 magmoms.append(read_float(line)) 

505 line = fd.readline() 

506 return magmoms, total_magmom 

507 

508 

509def read_energy(line, fd, debug=None): 

510 # It has Hartree unit yet 

511 return read_float(line) 

512 

513 

514def read_energies(line, fd, debug=None): 

515 line = fd.readline() 

516 if '***' in line: 

517 point = 7 # Version 3.8 

518 else: 

519 point = 16 # Version 3.9 

520 for _ in range(point): 

521 fd.readline() 

522 line = fd.readline() 

523 energies = [] 

524 while not (line == '' or line.isspace()): 

525 energies.append(float(line.split()[2])) 

526 line = fd.readline() 

527 return energies 

528 

529 

530def read_eigenvalues(line, fd, debug=False): 

531 """ 

532 Read the Eigenvalues in the `.out` file and returns the eigenvalue 

533 First, it assumes system have two spins and start reading until it reaches 

534 the end('*****...'). 

535 

536 eigenvalues[spin][kpoint][nbands] 

537 

538 For symmetry reason, `.out` file prints the eigenvalues at the half of the 

539 K points. Thus, we have to fill up the rest of the half. 

540 However, if the calculation was conducted only on the gamma point, it will 

541 raise the 'gamma_flag' as true and it will returns the original samples. 

542 """ 

543 def prind(*line, end='\n'): 

544 if debug: 

545 print(*line, end=end) 

546 prind("Read eigenvalues output") 

547 current_line = fd.tell() 

548 fd.seek(0) # Seek for the kgrid information 

549 while line != '': 

550 line = fd.readline().lower() 

551 if 'scf.kgrid' in line: 

552 break 

553 fd.seek(current_line) # Retrun to the original position 

554 

555 kgrid = read_tuple_integer(line) 

556 

557 if kgrid != (): 

558 prind('Non-Gamma point calculation') 

559 prind('scf.Kgrid is %d, %d, %d' % kgrid) 

560 gamma_flag = False 

561 # fd.seek(f.tell()+57) 

562 else: 

563 prind('Gamma point calculation') 

564 gamma_flag = True 

565 line = fd.readline() 

566 line = fd.readline() 

567 

568 eigenvalues = [] 

569 eigenvalues.append([]) 

570 eigenvalues.append([]) # Assume two spins 

571 i = 0 

572 while True: 

573 # Go to eigenvalues line 

574 while line != '': 

575 line = fd.readline() 

576 prind(line) 

577 ll = line.split() 

578 if line.isspace(): 

579 continue 

580 elif len(ll) > 1: 

581 if ll[0] == '1': 

582 break 

583 elif "*****" in line: 

584 break 

585 

586 # Break if it reaches the end or next parameters 

587 if "*****" in line or line == '': 

588 break 

589 

590 # Check Number and format is valid 

591 try: 

592 # Suppose to be looks like 

593 # 1 -2.33424746491277 -2.33424746917880 

594 ll = line.split() 

595 # Check if it reaches the end of the file 

596 assert line != '' 

597 assert len(ll) == 3 

598 float(ll[1]) 

599 float(ll[2]) 

600 except (AssertionError, ValueError): 

601 raise ParseError("Cannot read eigenvalues") 

602 

603 # Read eigenvalues 

604 eigenvalues[0].append([]) 

605 eigenvalues[1].append([]) 

606 while not (line == '' or line.isspace()): 

607 eigenvalues[0][i].append(float(rn(line, 2))) 

608 eigenvalues[1][i].append(float(rn(line, 1))) 

609 line = fd.readline() 

610 prind(line, end='') 

611 i += 1 

612 prind(line) 

613 if gamma_flag: 

614 return np.asarray(eigenvalues) 

615 eigen_half = np.asarray(eigenvalues) 

616 prind(eigen_half) 

617 # Fill up the half 

618 spin, half_kpts, bands = eigen_half.shape 

619 even_odd = np.array(kgrid).prod() % 2 

620 eigen_values = np.zeros((spin, half_kpts * 2 - even_odd, bands)) 

621 for i in range(half_kpts): 

622 eigen_values[0, i] = eigen_half[0, i, :] 

623 eigen_values[1, i] = eigen_half[1, i, :] 

624 eigen_values[0, 2 * half_kpts - 1 - i - even_odd] = eigen_half[0, i, :] 

625 eigen_values[1, 2 * half_kpts - 1 - i - even_odd] = eigen_half[1, i, :] 

626 return eigen_values 

627 

628 

629def read_forces(line, fd, debug=None): 

630 # It has Hartree per Bohr unit yet 

631 forces = [] 

632 fd.readline() # Skip Empty line 

633 line = fd.readline() 

634 while 'coordinates.forces>' not in line: 

635 forces.append(read_tuple_float(line)) 

636 line = fd.readline() 

637 return np.array(forces) 

638 

639 

640def read_dipole(line, fd, debug=None): 

641 dipole = [] 

642 while 'Total' not in line: 

643 line = fd.readline() 

644 dipole.append(read_tuple_float(line)) 

645 return dipole 

646 

647 

648def read_scaled_positions(line, fd, debug=None): 

649 scaled_positions = [] 

650 fd.readline() # Skip Empty lines 

651 fd.readline() 

652 fd.readline() 

653 line = fd.readline() 

654 while not (line == '' or line.isspace()): # Detect empty line 

655 scaled_positions.append(read_tuple_float(line)) 

656 line = fd.readline() 

657 return scaled_positions 

658 

659 

660def read_chemical_potential(line, fd, debug=None): 

661 return read_float(line) 

662 

663 

664def get_parameters(out_data=None, log_data=None, restart_data=None, 

665 scfout_data=None, dat_data=None, band_data=None): 

666 """ 

667 From the given data sets, construct the dictionary 'parameters'. If data 

668 is in the paramerters, it will save it. 

669 """ 

670 from ase.calculators.openmx import parameters as param 

671 scaned_data = [dat_data, out_data, log_data, restart_data, scfout_data, 

672 band_data] 

673 openmx_keywords = [param.tuple_integer_keys, param.tuple_float_keys, 

674 param.tuple_bool_keys, param.integer_keys, 

675 param.float_keys, param.string_keys, param.bool_keys, 

676 param.list_int_keys, param.list_bool_keys, 

677 param.list_float_keys, param.matrix_keys] 

678 parameters = {} 

679 for scaned_datum in scaned_data: 

680 for scaned_key in scaned_datum.keys(): 

681 for openmx_keyword in openmx_keywords: 

682 if scaned_key in get_standard_key(openmx_keyword): 

683 parameters[scaned_key] = scaned_datum[scaned_key] 

684 continue 

685 translated_parameters = get_standard_parameters(parameters) 

686 parameters.update(translated_parameters) 

687 return {k: v for k, v in parameters.items() if v is not None} 

688 

689 

690def get_standard_key(key): 

691 """ 

692 Standard ASE parameter format is to USE unerbar(_) instead of dot(.). Also, 

693 It is recommended to use lower case alphabet letter. Not Upper. Thus, we 

694 change the key to standard key 

695 For example: 

696 'scf.XcType' -> 'scf_xctype' 

697 """ 

698 if isinstance(key, str): 

699 return key.lower().replace('.', '_') 

700 elif isinstance(key, list): 

701 return [k.lower().replace('.', '_') for k in key] 

702 else: 

703 return [k.lower().replace('.', '_') for k in key] 

704 

705 

706def get_standard_parameters(parameters): 

707 """ 

708 Translate the OpenMX parameters to standard ASE parameters. For example, 

709 

710 scf.XcType -> xc 

711 scf.maxIter -> maxiter 

712 scf.energycutoff -> energy_cutoff 

713 scf.Kgrid -> kpts 

714 scf.EigenvalueSolver -> eigensolver 

715 scf.SpinPolarization -> spinpol 

716 scf.criterion -> convergence 

717 scf.Electric.Field -> external 

718 scf.Mixing.Type -> mixer 

719 scf.system.charge -> charge 

720 

721 We followed GPAW schem. 

722 """ 

723 from ase.calculators.openmx import parameters as param 

724 from ase.units import Bohr, Ha, Ry, fs, m, s 

725 units = param.unit_dat_keywords 

726 standard_parameters = {} 

727 standard_units = {'eV': 1, 'Ha': Ha, 'Ry': Ry, 'Bohr': Bohr, 'fs': fs, 

728 'K': 1, 'GV / m': 1e9 / 1.6e-19 / m, 'Ha/Bohr': Ha / Bohr, 

729 'm/s': m / s, '_amu': 1, 'Tesla': 1} 

730 translated_parameters = { 

731 'scf.XcType': 'xc', 

732 'scf.maxIter': 'maxiter', 

733 'scf.energycutoff': 'energy_cutoff', 

734 'scf.Kgrid': 'kpts', 

735 'scf.EigenvalueSolver': 'eigensolver', 

736 'scf.SpinPolarization': 'spinpol', 

737 'scf.criterion': 'convergence', 

738 'scf.Electric.Field': 'external', 

739 'scf.Mixing.Type': 'mixer', 

740 'scf.system.charge': 'charge' 

741 } 

742 

743 for key in parameters.keys(): 

744 for openmx_key, standard_key in translated_parameters.items(): 

745 if key == get_standard_key(openmx_key): 

746 unit = standard_units.get(units.get(openmx_key), 1) 

747 standard_parameters[standard_key] = parameters[key] * unit 

748 standard_parameters['spinpol'] = parameters.get('scf_spinpolarization') 

749 return standard_parameters 

750 

751 

752def get_atomic_formula(out_data=None, log_data=None, restart_data=None, 

753 scfout_data=None, dat_data=None): 

754 """_formula'. 

755 OpenMX results gives following information. Since, we should pick one