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

1# fmt: off 

2 

3# flake8: noqa 

4""" 

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

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

7functional theories. 

8 Copyright (C) 2018 JaeHwan Shim and JaeJun Yu 

9 

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

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

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

13 (at your option) any later version. 

14 

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

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

17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

18 GNU Lesser General Public License for more details. 

19 

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

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

22""" 

23# from ase.calculators import SinglePointDFTCalculator 

24import os 

25import struct 

26 

27import numpy as np 

28 

29from ase.io import ParseError 

30from ase.units import Bohr, Debye, Ha 

31 

32 

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

34 from ase import Atoms 

35 from ase.calculators.openmx import OpenMX 

36 """ 

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

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

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

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

41 next calculation results. 

42 

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

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

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

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

47 """ 

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

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

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

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

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

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

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

55 """ 

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

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

58 Fix the unit to ASE format. 

59 """ 

60 parameters = get_parameters(out_data=out_data, log_data=log_data, 

61 restart_data=restart_data, dat_data=dat_data, 

62 scfout_data=scfout_data, band_data=band_data) 

63 atomic_formula = get_atomic_formula(out_data=out_data, log_data=log_data, 

64 restart_data=restart_data, 

65 scfout_data=scfout_data, 

66 dat_data=dat_data) 

67 results = get_results(out_data=out_data, log_data=log_data, 

68 restart_data=restart_data, scfout_data=scfout_data, 

69 dat_data=dat_data, band_data=band_data) 

70 

71 atoms = Atoms(**atomic_formula) 

72 

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

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

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

76 from ase.calculators.calculator import all_properties 

77 from ase.calculators.singlepoint import SinglePointDFTCalculator 

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

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

80 

81 atoms.calc = SinglePointDFTCalculator( 

82 atoms, **results) 

83 

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

85 # atoms.calc.results = results 

86 return atoms 

87 

88 

89def read_file(filename, debug=False): 

90 """ 

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

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

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

94 the pair value of that key. For example, 

95 example will be written later 

96 """ 

97 from ase.calculators.openmx import parameters as param 

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

99 return {} 

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

101 'list_int_keys', 'list_float_keys', 'list_bool_keys', 

102 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys'] 

103 patterns = { 

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

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

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

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

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

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

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

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

112 special_patterns = { 

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

114 read_magmoms_and_total_magmom), 

115 } 

116 out_data = {} 

117 line = '\n' 

118 if (debug): 

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

120 with open(filename) as fd: 

121 ''' 

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

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

124 

125 if line in param.string_keys: 

126 out_data[key] = read_string(line) 

127 

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

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

130 

131 for [dtype] in parameters_keys: 

132 if line in param.[dtype]_keys: 

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

134 

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

136 ''' 

137 while line != '': 

138 pattern_matched = False 

139 line = fd.readline() 

140 try: 

141 _line = line.split()[0] 

142 except IndexError: 

143 continue 

144 for dtype_key in param_keys: 

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

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

147 for key in param.__dict__[dtype_key]: 

148 if key in _line: 

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

150 pattern_matched = True 

151 continue 

152 if pattern_matched: 

153 continue 

154 

155 for key in param.matrix_keys: 

156 if '<' + key in line: 

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

158 pattern_matched = True 

159 continue 

160 if pattern_matched: 

161 continue 

162 for key in patterns: 

163 if key in line: 

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

165 line, fd, debug=debug) 

166 pattern_matched = True 

167 continue 

168 if pattern_matched: 

169 continue 

170 for key in special_patterns: 

171 if key in line: 

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

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

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

175 pattern_matched = True 

176 continue 

177 if pattern_matched: 

178 continue 

179 return out_data 

180 

181 

182def read_scfout_file(filename=None): 

183 """ 

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

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

186 1, not 0 

187 

188 atomnum: the number of total atoms 

189 Catomnum: the number of atoms in the central region 

190 Latomnum: the number of atoms in the left lead 

191 Ratomnum: the number of atoms in the left lead 

192 SpinP_switch: 

193 0: non-spin polarized 

194 1: spin polarized 

195 TCpyCell: the total number of periodic cells 

196 Solver: method for solving eigenvalue problem 

197 ChemP: chemical potential 

198 Valence_Electrons: total number of valence electrons 

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

200 E_Temp: electronic temperature 

201 Total_NumOrbs: the number of atomic orbitals in each atom 

202 size: Total_NumOrbs[atomnum+1] 

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

204 size: FNAN[atomnum+1] 

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

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

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

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

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

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

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

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

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

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

215 note: 

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

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

218 Hks: Kohn-Sham matrix elements of basis orbitals 

219 size: Hks[SpinP_switch+1] 

220 [atomnum+1] 

221 [FNAN[ct_AN]+1] 

222 [Total_NumOrbs[ct_AN]] 

223 [Total_NumOrbs[h_AN]] 

224 iHks: 

225 imaginary Kohn-Sham matrix elements of basis orbitals 

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

227 of which contributions come from spin-orbit coupling 

228 and Hubbard U effective potential. 

229 size: iHks[3] 

230 [atomnum+1] 

231 [FNAN[ct_AN]+1] 

232 [Total_NumOrbs[ct_AN]] 

233 [Total_NumOrbs[h_AN]] 

234 OLP: overlap matrix 

235 size: OLP[atomnum+1] 

236 [FNAN[ct_AN]+1] 

237 [Total_NumOrbs[ct_AN]] 

238 [Total_NumOrbs[h_AN]] 

239 OLPpox: overlap matrix with position operator x 

240 size: OLPpox[atomnum+1] 

241 [FNAN[ct_AN]+1] 

242 [Total_NumOrbs[ct_AN]] 

243 [Total_NumOrbs[h_AN]] 

244 OLPpoy: overlap matrix with position operator y 

245 size: OLPpoy[atomnum+1] 

246 [FNAN[ct_AN]+1] 

247 [Total_NumOrbs[ct_AN]] 

248 [Total_NumOrbs[h_AN]] 

249 OLPpoz: overlap matrix with position operator z 

250 size: OLPpoz[atomnum+1] 

251 [FNAN[ct_AN]+1] 

252 [Total_NumOrbs[ct_AN]] 

253 [Total_NumOrbs[h_AN]] 

254 DM: overlap matrix 

255 size: DM[SpinP_switch+1] 

256 [atomnum+1] 

257 [FNAN[ct_AN]+1] 

258 [Total_NumOrbs[ct_AN]] 

259 [Total_NumOrbs[h_AN]] 

260 dipole_moment_core[4]: 

261 dipole_moment_background[4]: 

262 """ 

263 from numpy import cumsum as cum 

264 from numpy import insert as ins 

265 from numpy import split as spl 

266 from numpy import sum, zeros 

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

268 return {} 

269 

270 def easyReader(byte, data_type, shape): 

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

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

273 dt = data_type 

274 ds = data_size[data_type] 

275 unpack = struct.unpack 

276 if len(byte) == ds: 

277 if dt == 'i': 

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

279 elif dt == 'd': 

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

281 elif shape is not None: 

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

283 else: 

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

285 

286 def inte(byte, shape=None): 

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

288 

289 def floa(byte, shape=None): 

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

291 

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

293 myOLP = [] 

294 myOLP.append([]) 

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

296 myOLP.append([]) 

297 TNO1 = Total_NumOrbs[ct_AN] 

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

299 myOLP[ct_AN].append([]) 

300 Gh_AN = natn[ct_AN][h_AN] 

301 TNO2 = Total_NumOrbs[Gh_AN] 

302 for i in range(TNO1): 

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

304 return myOLP 

305 

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

307 Hks = [] 

308 for spin in range(SpinP_switch + 1): 

309 Hks.append([]) 

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

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

312 Hks[spin].append([]) 

313 TNO1 = Total_NumOrbs[ct_AN] 

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

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

316 Gh_AN = natn[ct_AN][h_AN] 

317 TNO2 = Total_NumOrbs[Gh_AN] 

318 for i in range(TNO1): 

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

320 return Hks 

321 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

340 iHks = [] 

341 if SpinP_switch == 3: 

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

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

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

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

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

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

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

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

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

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

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

353 

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

355 'Catomnum': Catomnum, 'Latomnum': Latomnum, 'Hks': Hks, 

356 'Ratomnum': Ratomnum, 'TCpyCell': TCpyCell, 'atv': atv, 

357 'Total_NumOrbs': Total_NumOrbs, 'FNAN': FNAN, 'natn': natn, 

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

359 'OLPpox': OLPpox, 'OLPpoy': OLPpoy, 'OLPpoz': OLPpoz, 

360 'Solver': Solver, 'ChemP': ChemP, 'E_Temp': E_Temp, 

361 'dipole_moment_core': dipole_moment_core, 'iHks': iHks, 

362 'dipole_moment_background': dipole_moment_background, 

363 'Valence_Electrons': Valence_Electrons, 'atv_ijk': atv_ijk, 

364 'Total_SpinS': Total_SpinS, 'DM': DM 

365 } 

366 return scf_out 

367 

368 

369def read_band_file(filename=None): 

370 band_data = {} 

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

372 return {} 

373 band_kpath = [] 

374 eigen_bands = [] 

375 with open(filename) as fd: 

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

377 nkpts = 0 

378 nband = int(line[0]) 

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

380 band_data['nband'] = nband 

381 band_data['nspin'] = nspin 

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

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

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

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

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

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

388 band_kpath.append(line) 

389 nkpts += int(line[0]) 

390 band_data['nkpts'] = nkpts 

391 band_data['band_kpath'] = band_kpath 

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

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

394 for i in range(nspin): 

395 for j in range(nkpts): 

396 line = fd.readline() 

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

398 line = fd.readline() 

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

400 band_data['eigenvalues'] = eigen_bands 

401 band_data['band_kpts'] = kpts 

402 return band_data 

403 

404 

405def read_electron_valency(filename='H_CA13'): 

406 array = [] 

407 with open(filename) as fd: 

408 array = fd.readlines() 

409 fd.close() 

410 required_line = '' 

411 for line in array: 

412 if 'valence.electron' in line: 

413 required_line = line 

414 return rn(required_line) 

415 

416 

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

418 """ 

419 Read n'th to last value. 

420 For example: 

421 ... 

422 scf.XcType LDA 

423 scf.Kgrid 4 4 4 

424 ... 

425 In Python, 

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

427 LDA 

428 >>> line = fd.readline() 

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

430 4 

431 """ 

432 return line.split()[-n] 

433 

434 

435def read_tuple_integer(line): 

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

437 

438 

439def read_tuple_float(line): 

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

441 

442 

443def read_integer(line): 

444 return int(rn(line)) 

445 

446 

447def read_float(line): 

448 return float(rn(line)) 

449 

450 

451def read_string(line): 

452 return str(rn(line)) 

453 

454 

455def read_bool(line): 

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

457 if bool == 'on': 

458 return True 

459 elif bool == 'off': 

460 return False 

461 else: 

462 return None 

463 

464 

465def read_list_int(line): 

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

467 

468 

469def read_list_float(line): 

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

471 

472 

473def read_list_bool(line): 

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

475 

476 

477def read_matrix(line, key, fd): 

478 matrix = [] 

479 line = fd.readline() 

480 while key not in line: 

481 matrix.append(line.split()) 

482 line = fd.readline() 

483 return matrix 

484 

485 

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

487 fd.readline() # passing empty line 

488 fd.readline() 

489 line = fd.readline() 

490 xx, xy, xz = read_tuple_float(line) 

491 line = fd.readline() 

492 yx, yy, yz = read_tuple_float(line) 

493 line = fd.readline() 

494 zx, zy, zz = read_tuple_float(line) 

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

496 return stress 

497 

498 

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

500 total_magmom = read_float(line) 

501 fd.readline() # Skip empty lines 

502 fd.readline() 

503 line = fd.readline() 

504 magmoms = [] 

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

506 magmoms.append(read_float(line)) 

507 line = fd.readline() 

508 return magmoms, total_magmom 

509 

510 

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

512 # It has Hartree unit yet 

513 return read_float(line) 

514 

515 

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

517 line = fd.readline() 

518 if '***' in line: 

519 point = 7 # Version 3.8 

520 else: 

521 point = 16 # Version 3.9 

522 for _ in range(point): 

523 fd.readline() 

524 line = fd.readline() 

525 energies = [] 

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

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

528 line = fd.readline() 

529 return energies 

530 

531 

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

533 """ 

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

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

536 the end('*****...'). 

537 

538 eigenvalues[spin][kpoint][nbands] 

539 

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

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

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

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

544 """ 

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

546 if debug: 

547 print(*line, end=end) 

548 prind("Read eigenvalues output") 

549 current_line = fd.tell() 

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

551 while line != '': 

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

553 if 'scf.kgrid' in line: 

554 break 

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

556 

557 kgrid = read_tuple_integer(line) 

558 

559 if kgrid != (): 

560 prind('Non-Gamma point calculation') 

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

562 gamma_flag = False 

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

564 else: 

565 prind('Gamma point calculation') 

566 gamma_flag = True 

567 line = fd.readline() 

568 line = fd.readline() 

569 

570 eigenvalues = [] 

571 eigenvalues.append([]) 

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

573 i = 0 

574 while True: 

575 # Go to eigenvalues line 

576 while line != '': 

577 line = fd.readline() 

578 prind(line) 

579 ll = line.split() 

580 if line.isspace(): 

581 continue 

582 elif len(ll) > 1: 

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

584 break 

585 elif "*****" in line: 

586 break 

587 

588 # Break if it reaches the end or next parameters 

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

590 break 

591 

592 # Check Number and format is valid 

593 try: 

594 # Suppose to be looks like 

595 # 1 -2.33424746491277 -2.33424746917880 

596 ll = line.split() 

597 # Check if it reaches the end of the file 

598 assert line != '' 

599 assert len(ll) == 3 

600 float(ll[1]) 

601 float(ll[2]) 

602 except (AssertionError, ValueError): 

603 raise ParseError("Cannot read eigenvalues") 

604 

605 # Read eigenvalues 

606 eigenvalues[0].append([]) 

607 eigenvalues[1].append([]) 

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

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

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

611 line = fd.readline() 

612 prind(line, end='') 

613 i += 1 

614 prind(line) 

615 if gamma_flag: 

616 return np.asarray(eigenvalues) 

617 eigen_half = np.asarray(eigenvalues) 

618 prind(eigen_half) 

619 # Fill up the half 

620 spin, half_kpts, bands = eigen_half.shape 

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

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

623 for i in range(half_kpts): 

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

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

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

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

628 return eigen_values 

629 

630 

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

632 # It has Hartree per Bohr unit yet 

633 forces = [] 

634 fd.readline() # Skip Empty line 

635 line = fd.readline() 

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

637 forces.append(read_tuple_float(line)) 

638 line = fd.readline() 

639 return np.array(forces) 

640 

641 

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

643 dipole = [] 

644 while 'Total' not in line: 

645 line = fd.readline() 

646 dipole.append(read_tuple_float(line)) 

647 return dipole 

648 

649 

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

651 scaled_positions = [] 

652 fd.readline() # Skip Empty lines 

653 fd.readline() 

654 fd.readline() 

655 line = fd.readline() 

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

657 scaled_positions.append(read_tuple_float(line)) 

658 line = fd.readline() 

659 return scaled_positions 

660 

661 

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

663 return read_float(line) 

664 

665 

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

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

668 """ 

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

670 is in the paramerters, it will save it. 

671 """ 

672 from ase.calculators.openmx import parameters as param 

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

674 band_data] 

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

676 param.tuple_bool_keys, param.integer_keys, 

677 param.float_keys, param.string_keys, param.bool_keys, 

678 param.list_int_keys, param.list_bool_keys, 

679 param.list_float_keys, param.matrix_keys] 

680 parameters = {} 

681 for scaned_datum in scaned_data: 

682 for scaned_key in scaned_datum.keys(): 

683 for openmx_keyword in openmx_keywords: 

684 if scaned_key in get_standard_key(openmx_keyword): 

685 parameters[scaned_key] = scaned_datum[scaned_key] 

686 continue 

687 translated_parameters = get_standard_parameters(parameters) 

688 parameters.update(translated_parameters) 

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

690 

691 

692def get_standard_key(key): 

693 """ 

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

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

696 change the key to standard key 

697 For example: 

698 'scf.XcType' -> 'scf_xctype' 

699 """ 

700 if isinstance(key, str): 

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

702 elif isinstance(key, list): 

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

704 else: 

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

706 

707 

708def get_standard_parameters(parameters): 

709 """ 

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

711 

712 scf.XcType -> xc 

713 scf.maxIter -> maxiter 

714 scf.energycutoff -> energy_cutoff 

715 scf.Kgrid -> kpts 

716 scf.EigenvalueSolver -> eigensolver 

717 scf.SpinPolarization -> spinpol 

718 scf.criterion -> convergence 

719 scf.Electric.Field -> external 

720 scf.Mixing.Type -> mixer 

721 scf.system.charge -> charge 

722 

723 We followed GPAW schem. 

724 """ 

725 from ase.calculators.openmx import parameters as param 

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

727 units = param.unit_dat_keywords 

728 standard_parameters = {} 

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

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

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

732 translated_parameters = { 

733 'scf.XcType': 'xc', 

734 'scf.maxIter': 'maxiter', 

735 'scf.energycutoff': 'energy_cutoff', 

736 'scf.Kgrid': 'kpts', 

737 'scf.EigenvalueSolver': 'eigensolver', 

738 'scf.SpinPolarization': 'spinpol', 

739 'scf.criterion': 'convergence', 

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

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

742 'scf.system.charge': 'charge' 

743 } 

744 

745 for key in parameters.keys(): 

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

747 if key == get_standard_key(openmx_key): 

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

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

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

751 return standard_parameters 

752 

753 

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

755 scfout_data=None, dat_data=None):