Coverage for /builds/kinetik161/ase/ase/calculators/openmx/reader.py: 63.28%

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 atoms.calc = OpenMX(**parameters) 

71 atoms.calc.results = results 

72 return atoms 

73 

74 

75def read_file(filename, debug=False): 

76 """ 

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

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

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

80 the pair value of that key. For example, 

81 example will be written later 

82 """ 

83 from ase.calculators.openmx import parameters as param 

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

85 return {} 

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

87 'list_int_keys', 'list_float_keys', 'list_bool_keys', 

88 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys'] 

89 patterns = { 

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

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

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

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

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

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

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

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

98 special_patterns = { 

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

100 read_magmoms_and_total_magmom), 

101 } 

102 out_data = {} 

103 line = '\n' 

104 if (debug): 

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

106 with open(filename) as fd: 

107 ''' 

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

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

110 

111 if line in param.string_keys: 

112 out_data[key] = read_string(line) 

113 

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

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

116 

117 for [dtype] in parameters_keys: 

118 if line in param.[dtype]_keys: 

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

120 

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

122 ''' 

123 while line != '': 

124 pattern_matched = False 

125 line = fd.readline() 

126 try: 

127 _line = line.split()[0] 

128 except IndexError: 

129 continue 

130 for dtype_key in param_keys: 

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

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

133 for key in param.__dict__[dtype_key]: 

134 if key in _line: 

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

136 pattern_matched = True 

137 continue 

138 if pattern_matched: 

139 continue 

140 

141 for key in param.matrix_keys: 

142 if '<' + key in line: 

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

144 pattern_matched = True 

145 continue 

146 if pattern_matched: 

147 continue 

148 for key in patterns: 

149 if key in line: 

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

151 line, fd, debug=debug) 

152 pattern_matched = True 

153 continue 

154 if pattern_matched: 

155 continue 

156 for key in special_patterns: 

157 if key in line: 

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

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

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

161 pattern_matched = True 

162 continue 

163 if pattern_matched: 

164 continue 

165 return out_data 

166 

167 

168def read_scfout_file(filename=None): 

169 """ 

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

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

172 1, not 0 

173 

174 atomnum: the number of total atoms 

175 Catomnum: the number of atoms in the central region 

176 Latomnum: the number of atoms in the left lead 

177 Ratomnum: the number of atoms in the left lead 

178 SpinP_switch: 

179 0: non-spin polarized 

180 1: spin polarized 

181 TCpyCell: the total number of periodic cells 

182 Solver: method for solving eigenvalue problem 

183 ChemP: chemical potential 

184 Valence_Electrons: total number of valence electrons 

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

186 E_Temp: electronic temperature 

187 Total_NumOrbs: the number of atomic orbitals in each atom 

188 size: Total_NumOrbs[atomnum+1] 

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

190 size: FNAN[atomnum+1] 

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

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

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

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

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

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

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

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

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

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

201 note: 

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

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

204 Hks: Kohn-Sham matrix elements of basis orbitals 

205 size: Hks[SpinP_switch+1] 

206 [atomnum+1] 

207 [FNAN[ct_AN]+1] 

208 [Total_NumOrbs[ct_AN]] 

209 [Total_NumOrbs[h_AN]] 

210 iHks: 

211 imaginary Kohn-Sham matrix elements of basis orbitals 

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

213 of which contributions come from spin-orbit coupling 

214 and Hubbard U effective potential. 

215 size: iHks[3] 

216 [atomnum+1] 

217 [FNAN[ct_AN]+1] 

218 [Total_NumOrbs[ct_AN]] 

219 [Total_NumOrbs[h_AN]] 

220 OLP: overlap matrix 

221 size: OLP[atomnum+1] 

222 [FNAN[ct_AN]+1] 

223 [Total_NumOrbs[ct_AN]] 

224 [Total_NumOrbs[h_AN]] 

225 OLPpox: overlap matrix with position operator x 

226 size: OLPpox[atomnum+1] 

227 [FNAN[ct_AN]+1] 

228 [Total_NumOrbs[ct_AN]] 

229 [Total_NumOrbs[h_AN]] 

230 OLPpoy: overlap matrix with position operator y 

231 size: OLPpoy[atomnum+1] 

232 [FNAN[ct_AN]+1] 

233 [Total_NumOrbs[ct_AN]] 

234 [Total_NumOrbs[h_AN]] 

235 OLPpoz: overlap matrix with position operator z 

236 size: OLPpoz[atomnum+1] 

237 [FNAN[ct_AN]+1] 

238 [Total_NumOrbs[ct_AN]] 

239 [Total_NumOrbs[h_AN]] 

240 DM: overlap matrix 

241 size: DM[SpinP_switch+1] 

242 [atomnum+1] 

243 [FNAN[ct_AN]+1] 

244 [Total_NumOrbs[ct_AN]] 

245 [Total_NumOrbs[h_AN]] 

246 dipole_moment_core[4]: 

247 dipole_moment_background[4]: 

248 """ 

249 from numpy import cumsum as cum 

250 from numpy import insert as ins 

251 from numpy import split as spl 

252 from numpy import sum, zeros 

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

254 return {} 

255 

256 def easyReader(byte, data_type, shape): 

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

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

259 dt = data_type 

260 ds = data_size[data_type] 

261 unpack = struct.unpack 

262 if len(byte) == ds: 

263 if dt == 'i': 

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

265 elif dt == 'd': 

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

267 elif shape is not None: 

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

269 else: 

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

271 

272 def inte(byte, shape=None): 

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

274 

275 def floa(byte, shape=None): 

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

277 

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

279 myOLP = [] 

280 myOLP.append([]) 

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

282 myOLP.append([]) 

283 TNO1 = Total_NumOrbs[ct_AN] 

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

285 myOLP[ct_AN].append([]) 

286 Gh_AN = natn[ct_AN][h_AN] 

287 TNO2 = Total_NumOrbs[Gh_AN] 

288 for i in range(TNO1): 

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

290 return myOLP 

291 

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

293 Hks = [] 

294 for spin in range(SpinP_switch + 1): 

295 Hks.append([]) 

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

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

298 Hks[spin].append([]) 

299 TNO1 = Total_NumOrbs[ct_AN] 

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

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

302 Gh_AN = natn[ct_AN][h_AN] 

303 TNO2 = Total_NumOrbs[Gh_AN] 

304 for i in range(TNO1): 

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

306 return Hks 

307 

308 fd = open(filename, mode='rb') 

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

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

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

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

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

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

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

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

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

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

319 tv = ins(floa(fd.read(8 * 3 * 4), shape=(3, 4)), 0, [0, 0, 0, 0], axis=0) 

320 rtv = ins(floa(fd.read(8 * 3 * 4), shape=(3, 4)), 0, [0, 0, 0, 0], axis=0) 

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

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

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

324 iHks = [] 

325 if SpinP_switch == 3: 

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

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

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

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

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

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

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

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

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

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

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

337 

338 fd.close() 

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

340 'Catomnum': Catomnum, 'Latomnum': Latomnum, 'Hks': Hks, 

341 'Ratomnum': Ratomnum, 'TCpyCell': TCpyCell, 'atv': atv, 

342 'Total_NumOrbs': Total_NumOrbs, 'FNAN': FNAN, 'natn': natn, 

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

344 'OLPpox': OLPpox, 'OLPpoy': OLPpoy, 'OLPpoz': OLPpoz, 

345 'Solver': Solver, 'ChemP': ChemP, 'E_Temp': E_Temp, 

346 'dipole_moment_core': dipole_moment_core, 'iHks': iHks, 

347 'dipole_moment_background': dipole_moment_background, 

348 'Valence_Electrons': Valence_Electrons, 'atv_ijk': atv_ijk, 

349 'Total_SpinS': Total_SpinS, 'DM': DM 

350 } 

351 return scf_out 

352 

353 

354def read_band_file(filename=None): 

355 band_data = {} 

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

357 return {} 

358 band_kpath = [] 

359 eigen_bands = [] 

360 with open(filename) as fd: 

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

362 nkpts = 0 

363 nband = int(line[0]) 

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

365 band_data['nband'] = nband 

366 band_data['nspin'] = nspin 

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

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

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

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

371 for i in range(band_data['band_nkpath']): 

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

373 band_kpath.append(line) 

374 nkpts += int(line[0]) 

375 band_data['nkpts'] = nkpts 

376 band_data['band_kpath'] = band_kpath 

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

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

379 for i in range(nspin): 

380 for j in range(nkpts): 

381 line = fd.readline() 

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

383 line = fd.readline() 

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

385 band_data['eigenvalues'] = eigen_bands 

386 band_data['band_kpts'] = kpts 

387 return band_data 

388 

389 

390def read_electron_valency(filename='H_CA13'): 

391 array = [] 

392 with open(filename) as fd: 

393 array = fd.readlines() 

394 fd.close() 

395 required_line = '' 

396 for line in array: 

397 if 'valence.electron' in line: 

398 required_line = line 

399 return rn(required_line) 

400 

401 

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

403 """ 

404 Read n'th to last value. 

405 For example: 

406 ... 

407 scf.XcType LDA 

408 scf.Kgrid 4 4 4 

409 ... 

410 In Python, 

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

412 LDA 

413 >>> line = fd.readline() 

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

415 4 

416 """ 

417 return line.split()[-n] 

418 

419 

420def read_tuple_integer(line): 

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

422 

423 

424def read_tuple_float(line): 

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

426 

427 

428def read_integer(line): 

429 return int(rn(line)) 

430 

431 

432def read_float(line): 

433 return float(rn(line)) 

434 

435 

436def read_string(line): 

437 return str(rn(line)) 

438 

439 

440def read_bool(line): 

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

442 if bool == 'on': 

443 return True 

444 elif bool == 'off': 

445 return False 

446 else: 

447 return None 

448 

449 

450def read_list_int(line): 

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

452 

453 

454def read_list_float(line): 

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

456 

457 

458def read_list_bool(line): 

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

460 

461 

462def read_matrix(line, key, fd): 

463 matrix = [] 

464 line = fd.readline() 

465 while key not in line: 

466 matrix.append(line.split()) 

467 line = fd.readline() 

468 return matrix 

469 

470 

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

472 fd.readline() # passing empty line 

473 fd.readline() 

474 line = fd.readline() 

475 xx, xy, xz = read_tuple_float(line) 

476 line = fd.readline() 

477 yx, yy, yz = read_tuple_float(line) 

478 line = fd.readline() 

479 zx, zy, zz = read_tuple_float(line) 

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

481 return stress 

482 

483 

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

485 total_magmom = read_float(line) 

486 fd.readline() # Skip empty lines 

487 fd.readline() 

488 line = fd.readline() 

489 magmoms = [] 

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

491 magmoms.append(read_float(line)) 

492 line = fd.readline() 

493 return magmoms, total_magmom 

494 

495 

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

497 # It has Hartree unit yet 

498 return read_float(line) 

499 

500 

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

502 line = fd.readline() 

503 if '***' in line: 

504 point = 7 # Version 3.8 

505 else: 

506 point = 16 # Version 3.9 

507 for i in range(point): 

508 fd.readline() 

509 line = fd.readline() 

510 energies = [] 

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

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

513 line = fd.readline() 

514 return energies 

515 

516 

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

518 """ 

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

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

521 the end('*****...'). 

522 

523 eigenvalues[spin][kpoint][nbands] 

524 

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

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

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

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

529 """ 

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

531 if debug: 

532 print(*line, end=end) 

533 prind("Read eigenvalues output") 

534 current_line = fd.tell() 

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

536 while line != '': 

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

538 if 'scf.kgrid' in line: 

539 break 

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

541 

542 kgrid = read_tuple_integer(line) 

543 

544 if kgrid != (): 

545 prind('Non-Gamma point calculation') 

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

547 gamma_flag = False 

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

549 else: 

550 prind('Gamma point calculation') 

551 gamma_flag = True 

552 line = fd.readline() 

553 line = fd.readline() 

554 

555 eigenvalues = [] 

556 eigenvalues.append([]) 

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

558 i = 0 

559 while True: 

560 # Go to eigenvalues line 

561 while line != '': 

562 line = fd.readline() 

563 prind(line) 

564 ll = line.split() 

565 if line.isspace(): 

566 continue 

567 elif len(ll) > 1: 

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

569 break 

570 elif "*****" in line: 

571 break 

572 

573 # Break if it reaches the end or next parameters 

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

575 break 

576 

577 # Check Number and format is valid 

578 try: 

579 # Suppose to be looks like 

580 # 1 -2.33424746491277 -2.33424746917880 

581 ll = line.split() 

582 # Check if it reaches the end of the file 

583 assert line != '' 

584 assert len(ll) == 3 

585 float(ll[1]) 

586 float(ll[2]) 

587 except (AssertionError, ValueError): 

588 raise ParseError("Cannot read eigenvalues") 

589 

590 # Read eigenvalues 

591 eigenvalues[0].append([]) 

592 eigenvalues[1].append([]) 

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

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

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

596 line = fd.readline() 

597 prind(line, end='') 

598 i += 1 

599 prind(line) 

600 if gamma_flag: 

601 return np.asarray(eigenvalues) 

602 eigen_half = np.asarray(eigenvalues) 

603 prind(eigen_half) 

604 # Fill up the half 

605 spin, half_kpts, bands = eigen_half.shape 

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

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

608 for i in range(half_kpts): 

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

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

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

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

613 return eigen_values 

614 

615 

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

617 # It has Hartree per Bohr unit yet 

618 forces = [] 

619 fd.readline() # Skip Empty line 

620 line = fd.readline() 

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

622 forces.append(read_tuple_float(line)) 

623 line = fd.readline() 

624 return np.array(forces) 

625 

626 

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

628 dipole = [] 

629 while 'Total' not in line: 

630 line = fd.readline() 

631 dipole.append(read_tuple_float(line)) 

632 return dipole 

633 

634 

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

636 scaled_positions = [] 

637 fd.readline() # Skip Empty lines 

638 fd.readline() 

639 fd.readline() 

640 line = fd.readline() 

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

642 scaled_positions.append(read_tuple_float(line)) 

643 line = fd.readline() 

644 return scaled_positions 

645 

646 

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

648 return read_float(line) 

649 

650 

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

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

653 """ 

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

655 is in the paramerters, it will save it. 

656 """ 

657 from ase.calculators.openmx import parameters as param 

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

659 band_data] 

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

661 param.tuple_bool_keys, param.integer_keys, 

662 param.float_keys, param.string_keys, param.bool_keys, 

663 param.list_int_keys, param.list_bool_keys, 

664 param.list_float_keys, param.matrix_keys] 

665 parameters = {} 

666 for scaned_datum in scaned_data: 

667 for scaned_key in scaned_datum.keys(): 

668 for openmx_keyword in openmx_keywords: 

669 if scaned_key in get_standard_key(openmx_keyword): 

670 parameters[scaned_key] = scaned_datum[scaned_key] 

671 continue 

672 translated_parameters = get_standard_parameters(parameters) 

673 parameters.update(translated_parameters) 

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

675 

676 

677def get_standard_key(key): 

678 """ 

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

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

681 change the key to standard key 

682 For example: 

683 'scf.XcType' -> 'scf_xctype' 

684 """ 

685 if isinstance(key, str): 

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

687 elif isinstance(key, list): 

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

689 else: 

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

691 

692 

693def get_standard_parameters(parameters): 

694 """ 

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

696 

697 scf.XcType -> xc 

698 scf.maxIter -> maxiter 

699 scf.energycutoff -> energy_cutoff 

700 scf.Kgrid -> kpts 

701 scf.EigenvalueSolver -> eigensolver 

702 scf.SpinPolarization -> spinpol 

703 scf.criterion -> convergence 

704 scf.Electric.Field -> external 

705 scf.Mixing.Type -> mixer 

706 scf.system.charge -> charge 

707 

708 We followed GPAW schem. 

709 """ 

710 from ase.calculators.openmx import parameters as param 

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

712 units = param.unit_dat_keywords 

713 standard_parameters = {} 

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

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

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

717 translated_parameters = { 

718 'scf.XcType': 'xc', 

719 'scf.maxIter': 'maxiter', 

720 'scf.energycutoff': 'energy_cutoff', 

721 'scf.Kgrid': 'kpts', 

722 'scf.EigenvalueSolver': 'eigensolver', 

723 'scf.SpinPolarization': 'spinpol', 

724 'scf.criterion': 'convergence', 

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

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

727 'scf.system.charge': 'charge' 

728 } 

729 

730 for key in parameters.keys(): 

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

732 if key == get_standard_key(openmx_key): 

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

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

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

736 return standard_parameters 

737 

738 

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

740 scfout_data=None, dat_data=None): 

741 """_formula'. 

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

743 between position/scaled_position, scaled_positions are suppressed by 

744 default. We use input value of position. Not the position after 

745 calculation. It is temporal. 

746 

747 Atoms.SpeciesAndCoordinate -> symbols 

748 Atoms.SpeciesAndCoordinate -> positions 

749 Atoms.UnitVectors -> cell 

750 scaled_positions -> scaled_positions 

751 If `positions` and `scaled_positions` are both given, this key deleted 

752 magmoms -> magmoms, Single value for each atom or three numbers for each 

753 atom for non-collinear calculations. 

754 """ 

755 atomic_formula = {} 

756 parameters = {'symbols': list, 'positions': list, 'scaled_positions': list,