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 24import numpy as np 25from ase.units import Ha, Bohr, Debye 26from ase.io import ParseError 27 28 29 30def read_openmx(filename=None, debug=False): 31 from ase.calculators.openmx import OpenMX 32 from ase import Atoms 33 """ 34 Read results from typical OpenMX output files and returns the atom object 35 In default mode, it reads every implementd properties we could get from 36 the files. Unlike previous version, we read the information based on file. 37 previous results will be eraised unless previous results are written in the 38 next calculation results. 39 40 Read the 'LABEL.log' file seems redundant. Because all the 41 information should already be written in '.out' file. However, in the 42 version 3.8.3, stress tensor are not written in the '.out' file. It only 43 contained in the '.log' file. So... I implented reading '.log' file method 44 """ 45 log_data = read_file(get_file_name('.log', filename), debug=debug) 46 restart_data = read_file(get_file_name('.dat#', filename), debug=debug) 47 dat_data = read_file(get_file_name('.dat', filename), debug=debug) 48 out_data = read_file(get_file_name('.out', filename), debug=debug) 49 scfout_data = read_scfout_file(get_file_name('.scfout', filename)) 50 band_data = read_band_file(get_file_name('.Band', filename)) 51 # dos_data = read_dos_file(get_file_name('.Dos.val', filename)) 52 """ 53 First, we get every data we could get from the all results files. And then, 54 reform the data to fit to data structure of Atom object. While doing this, 55 Fix the unit to ASE format. 56 """ 57 parameters = get_parameters(out_data=out_data, log_data=log_data, 58 restart_data=restart_data, dat_data=dat_data, 59 scfout_data=scfout_data, band_data=band_data) 60 atomic_formula = get_atomic_formula(out_data=out_data, log_data=log_data, 61 restart_data=restart_data, 62 scfout_data=scfout_data, 63 dat_data=dat_data) 64 results = get_results(out_data=out_data, log_data=log_data, 65 restart_data=restart_data, scfout_data=scfout_data, 66 dat_data=dat_data, band_data=band_data) 67 68 atoms = Atoms(**atomic_formula) 69 atoms.calc = OpenMX(**parameters) 70 atoms.calc.results = results 71 return atoms 72 73 74def read_file(filename, debug=False): 75 """ 76 Read the 'LABEL.out' file. Using 'parameters.py', we read every 'allowed_ 77 dat' dictionory. while reading a file, if one find the key matcheds That 78 'patters', which indicates the property we want is written, it will returns 79 the pair value of that key. For example, 80 example will be written later 81 """ 82 from ase.calculators.openmx import parameters as param 83 if not os.path.isfile(filename): 84 return {} 85 param_keys = ['integer_keys', 'float_keys', 'string_keys', 'bool_keys', 86 'list_int_keys', 'list_float_keys', 'list_bool_keys', 87 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys'] 88 patterns = { 89 'Stress tensor': ('stress', read_stress_tensor), 90 'Dipole moment': ('dipole', read_dipole), 91 'Fractional coordinates of': ('scaled_positions', read_scaled_positions), 92 'Utot.': ('energy', read_energy), 93 'energies in': ('energies', read_energies), 94 'Chemical Potential': ('chemical_potential', read_chemical_potential), 95 '<coordinates.forces': ('forces', read_forces), 96 'Eigenvalues (Hartree)': ('eigenvalues', read_eigenvalues)} 97 special_patterns = { 98 'Total spin moment': (('magmoms', 'total_magmom'), 99 read_magmoms_and_total_magmom), 100 } 101 out_data = {} 102 line = '\n' 103 if(debug): 104 print('Read results from %s' % filename) 105 with open(filename, 'r') as fd: 106 ''' 107 Read output file line by line. When the `line` matches the pattern 108 of certain keywords in `param.[dtype]_keys`, for example, 109 110 if line in param.string_keys: 111 out_data[key] = read_string(line) 112 113 parse that line and store it to `out_data` in specified data type. 114 To cover all `dtype` parameters, for loop was used, 115 116 for [dtype] in parameters_keys: 117 if line in param.[dtype]_keys: 118 out_data[key] = read_[dtype](line) 119 120 After found matched pattern, escape the for loop using `continue`. 121 ''' 122 while line != '': 123 pattern_matched = False 124 line = fd.readline() 125 try: 126 _line = line.split()[0] 127 except IndexError: 128 continue 129 for dtype_key in param_keys: 130 dtype = dtype_key.rsplit('_', 1)[0] 131 read_dtype = globals()['read_' + dtype] 132 for key in param.__dict__[dtype_key]: 133 if key in _line: 134 out_data[get_standard_key(key)] = read_dtype(line) 135 pattern_matched = True 136 continue 137 if pattern_matched: 138 continue 139 140 for key in param.matrix_keys: 141 if '<'+key in line: 142 out_data[get_standard_key(key)] = read_matrix(line, key, fd) 143 pattern_matched = True 144 continue 145 if pattern_matched: 146 continue 147 for key in patterns.keys(): 148 if key in line: 149 out_data[patterns[key][0]] = patterns[key][1](line, fd, debug=debug) 150 pattern_matched = True 151 continue 152 if pattern_matched: 153 continue 154 for key in special_patterns.keys(): 155 if key in line: 156 a, b = special_patterns[key][1](line, fd) 157 out_data[special_patterns[key][0][0]] = a 158 out_data[special_patterns[key][0][1]] = b 159 pattern_matched = True 160 continue 161 if pattern_matched: 162 continue 163 return out_data 164 165 166def read_scfout_file(filename=None): 167 """ 168 Read the Developer output '.scfout' files. It Behaves like read_scfout.c, 169 OpenMX module, but written in python. Note that some array are begin with 170 1, not 0 171 172 atomnum: the number of total atoms 173 Catomnum: the number of atoms in the central region 174 Latomnum: the number of atoms in the left lead 175 Ratomnum: the number of atoms in the left lead 176 SpinP_switch: 177 0: non-spin polarized 178 1: spin polarized 179 TCpyCell: the total number of periodic cells 180 Solver: method for solving eigenvalue problem 181 ChemP: chemical potential 182 Valence_Electrons: total number of valence electrons 183 Total_SpinS: total value of Spin (2*Total_SpinS = muB) 184 E_Temp: electronic temperature 185 Total_NumOrbs: the number of atomic orbitals in each atom 186 size: Total_NumOrbs[atomnum+1] 187 FNAN: the number of first neighboring atoms of each atom 188 size: FNAN[atomnum+1] 189 natn: global index of neighboring atoms of an atom ct_AN 190 size: natn[atomnum+1][FNAN[ct_AN]+1] 191 ncn: global index for cell of neighboring atoms of an atom ct_AN 192 size: ncn[atomnum+1][FNAN[ct_AN]+1] 193 atv: x,y,and z-components of translation vector of periodically copied cell 194 size: atv[TCpyCell+1][4]: 195 atv_ijk: i,j,and j number of periodically copied cells 196 size: atv_ijk[TCpyCell+1][4]: 197 tv[4][4]: unit cell vectors in Bohr 198 rtv[4][4]: reciprocal unit cell vectors in Bohr^{-1} 199 note: 200 tv_i dot rtv_j = 2PI * Kronecker's delta_{ij} 201 Gxyz[atomnum+1][60]: atomic coordinates in Bohr 202 Hks: Kohn-Sham matrix elements of basis orbitals 203 size: Hks[SpinP_switch+1] 204 [atomnum+1] 205 [FNAN[ct_AN]+1] 206 [Total_NumOrbs[ct_AN]] 207 [Total_NumOrbs[h_AN]] 208 iHks: 209 imaginary Kohn-Sham matrix elements of basis orbitals 210 for alpha-alpha, beta-beta, and alpha-beta spin matrices 211 of which contributions come from spin-orbit coupling 212 and Hubbard U effective potential. 213 size: iHks[3] 214 [atomnum+1] 215 [FNAN[ct_AN]+1] 216 [Total_NumOrbs[ct_AN]] 217 [Total_NumOrbs[h_AN]] 218 OLP: overlap matrix 219 size: OLP[atomnum+1] 220 [FNAN[ct_AN]+1] 221 [Total_NumOrbs[ct_AN]] 222 [Total_NumOrbs[h_AN]] 223 OLPpox: overlap matrix with position operator x 224 size: OLPpox[atomnum+1] 225 [FNAN[ct_AN]+1] 226 [Total_NumOrbs[ct_AN]] 227 [Total_NumOrbs[h_AN]] 228 OLPpoy: overlap matrix with position operator y 229 size: OLPpoy[atomnum+1] 230 [FNAN[ct_AN]+1] 231 [Total_NumOrbs[ct_AN]] 232 [Total_NumOrbs[h_AN]] 233 OLPpoz: overlap matrix with position operator z 234 size: OLPpoz[atomnum+1] 235 [FNAN[ct_AN]+1] 236 [Total_NumOrbs[ct_AN]] 237 [Total_NumOrbs[h_AN]] 238 DM: overlap matrix 239 size: DM[SpinP_switch+1] 240 [atomnum+1] 241 [FNAN[ct_AN]+1] 242 [Total_NumOrbs[ct_AN]] 243 [Total_NumOrbs[h_AN]] 244 dipole_moment_core[4]: 245 dipole_moment_background[4]: 246 """ 247 from numpy import insert as ins 248 from numpy import cumsum as cum 249 from numpy import split as spl 250 from numpy import sum, zeros 251 if not os.path.isfile(filename): 252 return {} 253 254 def easyReader(byte, data_type, shape): 255 data_size = {'d': 8, 'i': 4} 256 data_struct = {'d': float, 'i': int} 257 dt = data_type 258 ds = data_size[data_type] 259 unpack = struct.unpack 260 if len(byte) == ds: 261 if dt == 'i': 262 return data_struct[dt].from_bytes(byte, byteorder='little') 263 elif dt == 'd': 264 return np.array(unpack(dt*(len(byte)//ds), byte))[0] 265 elif shape is not None: 266 return np.array(unpack(dt*(len(byte)//ds), byte)).reshape(shape) 267 else: 268 return np.array(unpack(dt*(len(byte)//ds), byte)) 269 270 def inte(byte, shape=None): 271 return easyReader(byte, 'i', shape) 272 273 def floa(byte, shape=None): 274 return easyReader(byte, 'd', shape) 275 276 def readOverlap(atomnum, Total_NumOrbs, FNAN, natn, f): 277 myOLP = [] 278 myOLP.append([]) 279 for ct_AN in range(1, atomnum + 1): 280 myOLP.append([]) 281 TNO1 = Total_NumOrbs[ct_AN] 282 for h_AN in range(FNAN[ct_AN] + 1): 283 myOLP[ct_AN].append([]) 284 Gh_AN = natn[ct_AN][h_AN] 285 TNO2 = Total_NumOrbs[Gh_AN] 286 for i in range(TNO1): 287 myOLP[ct_AN][h_AN].append(floa(f.read(8*TNO2))) 288 return myOLP 289 290 def readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, f): 291 Hks = [] 292 for spin in range(SpinP_switch + 1): 293 Hks.append([]) 294 Hks[spin].append([np.zeros(FNAN[0] + 1)]) 295 for ct_AN in range(1, atomnum + 1): 296 Hks[spin].append([]) 297 TNO1 = Total_NumOrbs[ct_AN] 298 for h_AN in range(FNAN[ct_AN] + 1): 299 Hks[spin][ct_AN].append([]) 300 Gh_AN = natn[ct_AN][h_AN] 301 TNO2 = Total_NumOrbs[Gh_AN] 302 for i in range(TNO1): 303 Hks[spin][ct_AN][h_AN].append(floa(f.read(8*TNO2))) 304 return Hks 305 306 fd = open(filename, mode='rb') 307 atomnum, SpinP_switch = inte(fd.read(8)) 308 Catomnum, Latomnum, Ratomnum, TCpyCell = inte(fd.read(16)) 309 atv = floa(fd.read(8*4*(TCpyCell+1)), shape=(TCpyCell+1, 4)) 310 atv_ijk = inte(fd.read(4*4*(TCpyCell+1)), shape=(TCpyCell+1, 4)) 311 Total_NumOrbs = np.insert(inte(fd.read(4*(atomnum))), 0, 1, axis=0) 312 FNAN = np.insert(inte(fd.read(4*(atomnum))), 0, 0, axis=0) 313 natn = ins(spl(inte(fd.read(4*sum(FNAN[1:] + 1))), cum(FNAN[1:] + 1)), 314 0, zeros(FNAN[0] + 1), axis=0)[:-1] 315 ncn = ins(spl(inte(fd.read(4*np.sum(FNAN[1:] + 1))), cum(FNAN[1:] + 1)), 316 0, np.zeros(FNAN[0] + 1), axis=0)[:-1] 317 tv = ins(floa(fd.read(8*3*4), shape=(3, 4)), 0, [0, 0, 0, 0], axis=0) 318 rtv = ins(floa(fd.read(8*3*4), shape=(3, 4)), 0, [0, 0, 0, 0], axis=0) 319 Gxyz = ins(floa(fd.read(8*(atomnum)*4), shape=(atomnum, 4)), 0, 320 [0., 0., 0., 0.], axis=0) 321 Hks = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd) 322 iHks = [] 323 if SpinP_switch == 3: 324 iHks = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd) 325 OLP = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 326 OLPpox = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 327 OLPpoy = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 328 OLPpoz = readOverlap(atomnum, Total_NumOrbs, FNAN, natn, fd) 329 DM = readHam(SpinP_switch, FNAN, atomnum, Total_NumOrbs, natn, fd) 330 Solver = inte(fd.read(4)) 331 ChemP, E_Temp = floa(fd.read(8*2)) 332 dipole_moment_core = floa(fd.read(8*3)) 333 dipole_moment_background = floa(fd.read(8*3)) 334 Valence_Electrons, Total_SpinS = floa(fd.read(8*2)) 335 336 fd.close() 337 scf_out = {'atomnum': atomnum, 'SpinP_switch': SpinP_switch, 338 'Catomnum': Catomnum, 'Latomnum': Latomnum, 'Hks': Hks, 339 'Ratomnum': Ratomnum, 'TCpyCell': TCpyCell, 'atv': atv, 340 'Total_NumOrbs': Total_NumOrbs, 'FNAN': FNAN, 'natn': natn, 341 'ncn': ncn, 'tv': tv, 'rtv': rtv, 'Gxyz': Gxyz, 'OLP': OLP, 342 'OLPpox': OLPpox, 'OLPpoy': OLPpoy, 'OLPpoz': OLPpoz, 343 'Solver': Solver, 'ChemP': ChemP, 'E_Temp': E_Temp, 344 'dipole_moment_core': dipole_moment_core, 'iHks': iHks, 345 'dipole_moment_background': dipole_moment_background, 346 'Valence_Electrons': Valence_Electrons, 'atv_ijk': atv_ijk, 347 'Total_SpinS': Total_SpinS, 'DM': DM 348 } 349 return scf_out 350 351 352def read_band_file(filename=None): 353 band_data = {} 354 if not os.path.isfile(filename): 355 return {} 356 band_kpath = [] 357 eigen_bands = [] 358 with open(filename, 'r') as fd: 359 line = f.readline().split() 360 nkpts = 0 361 nband = int(line[0]) 362 nspin = int(line[1]) + 1 363 band_data['nband'] = nband 364 band_data['nspin'] = nspin 365 line = f.readline().split() 366 band_data['band_kpath_unitcell'] = [line[:3], line[3:6], line[6:9]] 367 line = f.readline().split() 368 band_data['band_nkpath'] = int(line[0]) 369 for i in range(band_data['band_nkpath']): 370 line = f.readline().split() 371 band_kpath.append(line) 372 nkpts += int(line[0]) 373 band_data['nkpts'] = nkpts 374 band_data['band_kpath'] = band_kpath 375 kpts = np.zeros((nkpts, 3)) 376 eigen_bands = np.zeros((nspin, nkpts, nband)) 377 for i in range(nspin): 378 for j in range(nkpts): 379 line = f.readline() 380 kpts[j] = np.array(line.split(), dtype=float)[1:] 381 line = f.readline() 382 eigen_bands[i, j] = np.array(line.split(), dtype=float)[:] 383 band_data['eigenvalues'] = eigen_bands 384 band_data['band_kpts'] = kpts 385 return band_data 386 387 388def read_electron_valency(filename='H_CA13'): 389 array = [] 390 with open(filename, 'r') as fd: 391 array = fd.readlines() 392 fd.close() 393 required_line = '' 394 for line in array: 395 if 'valence.electron' in line: 396 required_line = line 397 return rn(required_line) 398 399 400def rn(line='\n', n=1): 401 """ 402 Read n'th to last value. 403 For example: 404 ... 405 scf.XcType LDA 406 scf.Kgrid 4 4 4 407 ... 408 In Python, 409 >>> str(rn(line, 1)) 410 LDA 411 >>> line = f.readline() 412 >>> int(rn(line, 3)) 413 4 414 """ 415 return line.split()[-n] 416 417 418def read_tuple_integer(line): 419 return tuple([int(x) for x in line.split()[-3:]]) 420 421 422def read_tuple_float(line): 423 return tuple([float(x) for x in line.split()[-3:]]) 424 425 426def read_integer(line): 427 return int(rn(line)) 428 429 430def read_float(line): 431 return float(rn(line)) 432 433 434def read_string(line): 435 return str(rn(line)) 436 437 438def read_bool(line): 439 bool = str(rn(line)).lower() 440 if bool == 'on': 441 return True 442 elif bool == 'off': 443 return False 444 else: 445 return None 446 447 448def read_list_int(line): 449 return [int(x) for x in line.split()[1:]] 450 451 452def read_list_float(line): 453 return [float(x) for x in line.split()[1:]] 454 455 456def read_list_bool(line): 457 return [read_bool(x) for x in line.split()[1:]] 458 459 460def read_matrix(line, key, f): 461 matrix = [] 462 line = f.readline() 463 while key not in line: 464 matrix.append(line.split()) 465 line = f.readline() 466 return matrix 467 468 469def read_stress_tensor(line, f, debug=None): 470 f.readline() # passing empty line 471 f.readline() 472 line = f.readline() 473 xx, xy, xz = read_tuple_float(line) 474 line = f.readline() 475 yx, yy, yz = read_tuple_float(line) 476 line = f.readline() 477 zx, zy, zz = read_tuple_float(line) 478 stress = [xx, yy, zz, (zy + yz)/2, (zx + xz)/2, (yx + xy)/2] 479 return stress 480 481 482def read_magmoms_and_total_magmom(line, f, debug=None): 483 total_magmom = read_float(line) 484 f.readline() # Skip empty lines 485 f.readline() 486 line = f.readline() 487 magmoms = [] 488 while not(line == '' or line.isspace()): 489 magmoms.append(read_float(line)) 490 line = f.readline() 491 return magmoms, total_magmom 492 493 494def read_energy(line, f, debug=None): 495 # It has Hartree unit yet 496 return read_float(line) 497 498def read_energies(line, f, debug=None): 499 line = f.readline() 500 if '***' in line: 501 point = 7 # Version 3.8 502 else: 503 point = 16 # Version 3.9 504 for i in range(point): 505 f.readline() 506 line = f.readline() 507 energies = [] 508 while not(line == '' or line.isspace()): 509 energies.append(float(line.split()[2])) 510 line = f.readline() 511 return energies 512 513def read_eigenvalues(line, f, debug=False): 514 """ 515 Read the Eigenvalues in the `.out` file and returns the eigenvalue 516 First, it assumes system have two spins and start reading until it reaches 517 the end('*****...'). 518 519 eigenvalues[spin][kpoint][nbands] 520 521 For symmetry reason, `.out` file prints the eigenvalues at the half of the 522 K points. Thus, we have to fill up the rest of the half. 523 However, if the calculation was conducted only on the gamma point, it will 524 raise the 'gamma_flag' as true and it will returns the original samples. 525 """ 526 def prind(*line, end='\n'): 527 if debug: 528 print(*line, end=end) 529 prind("Read eigenvalues output") 530 current_line = f.tell() 531 f.seek(0) # Seek for the kgrid information 532 while line != '': 533 line = f.readline().lower() 534 if 'scf.kgrid' in line: 535 break 536 f.seek(current_line) # Retrun to the original position 537 538 kgrid = read_tuple_integer(line) 539 540 if kgrid != (): 541 prind('Non-Gamma point calculation') 542 prind('scf.Kgrid is %d, %d, %d' % kgrid) 543 gamma_flag = False 544 # f.seek(f.tell()+57) 545 else: 546 prind('Gamma point calculation') 547 gamma_flag = True 548 line = f.readline() 549 line = f.readline() 550 551 eigenvalues = [] 552 eigenvalues.append([]) 553 eigenvalues.append([]) # Assume two spins 554 i = 0 555 while True: 556 # Go to eigenvalues line 557 while line != '': 558 line = f.readline() 559 prind(line) 560 ll = line.split() 561 if line.isspace(): 562 continue 563 elif len(ll) > 1: 564 if ll[0] == '1': 565 break 566 elif "*****" in line: 567 break 568 569 # Break if it reaches the end or next parameters 570 if "*****" in line or line == '': 571 break 572 573 # Check Number and format is valid 574 try: 575 # Suppose to be looks like 576 # 1 -2.33424746491277 -2.33424746917880 577 ll = line.split() 578 # Check if it reaches the end of the file 579 assert line != '' 580 assert len(ll) == 3 581 float(ll[1]); float(ll[2]) 582 except (AssertionError, ValueError): 583 raise ParseError("Cannot read eigenvalues") 584 585 # Read eigenvalues 586 eigenvalues[0].append([]) 587 eigenvalues[1].append([]) 588 while not (line == '' or line.isspace()): 589 eigenvalues[0][i].append(float(rn(line, 2))) 590 eigenvalues[1][i].append(float(rn(line, 1))) 591 line = f.readline() 592 prind(line, end='') 593 i += 1 594 prind(line) 595 if gamma_flag: 596 return np.asarray(eigenvalues) 597 eigen_half = np.asarray(eigenvalues) 598 prind(eigen_half) 599 # Fill up the half 600 spin, half_kpts, bands = eigen_half.shape 601 even_odd = np.array(kgrid).prod() % 2 602 eigen_values = np.zeros((spin, half_kpts*2-even_odd, bands)) 603 for i in range(half_kpts): 604 eigen_values[0, i] = eigen_half[0, i, :] 605 eigen_values[1, i] = eigen_half[1, i, :] 606 eigen_values[0, 2*half_kpts-1-i-even_odd] = eigen_half[0, i, :] 607 eigen_values[1, 2*half_kpts-1-i-even_odd] = eigen_half[1, i, :] 608 return eigen_values 609 610 611def read_forces(line, f, debug=None): 612 # It has Hartree per Bohr unit yet 613 forces = [] 614 f.readline() # Skip Empty line 615 line = f.readline() 616 while 'coordinates.forces>' not in line: 617 forces.append(read_tuple_float(line)) 618 line = f.readline() 619 return np.array(forces) 620 621 622def read_dipole(line, f, debug=None): 623 dipole = [] 624 while 'Total' not in line: 625 line = f.readline() 626 dipole.append(read_tuple_float(line)) 627 return dipole 628 629 630def read_scaled_positions(line, f, debug=None): 631 scaled_positions = [] 632 f.readline() # Skip Empty lines 633 f.readline() 634 f.readline() 635 line = f.readline() 636 while not(line == '' or line.isspace()): # Detect empty line 637 scaled_positions.append(read_tuple_float(line)) 638 line = f.readline() 639 return scaled_positions 640 641 642def read_chemical_potential(line, f, debug=None): 643 return read_float(line) 644 645 646def get_parameters(out_data=None, log_data=None, restart_data=None, 647 scfout_data=None, dat_data=None, band_data=None): 648 """ 649 From the given data sets, construct the dictionary 'parameters'. If data 650 is in the paramerters, it will save it. 651 """ 652 from ase.calculators.openmx import parameters as param 653 scaned_data = [dat_data, out_data, log_data, restart_data, scfout_data, 654 band_data] 655 openmx_keywords = [param.tuple_integer_keys, param.tuple_float_keys, 656 param.tuple_bool_keys, param.integer_keys, 657 param.float_keys, param.string_keys, param.bool_keys, 658 param.list_int_keys, param.list_bool_keys, 659 param.list_float_keys, param.matrix_keys] 660 parameters = {} 661 for scaned_datum in scaned_data: 662 for scaned_key in scaned_datum.keys(): 663 for openmx_keyword in openmx_keywords: 664 if scaned_key in get_standard_key(openmx_keyword): 665 parameters[scaned_key] = scaned_datum[scaned_key] 666 continue 667 translated_parameters = get_standard_parameters(parameters) 668 parameters.update(translated_parameters) 669 return {k: v for k, v in parameters.items() if v is not None} 670 671 672def get_standard_key(key): 673 """ 674 Standard ASE parameter format is to USE unerbar(_) instead of dot(.). Also, 675 It is recommended to use lower case alphabet letter. Not Upper. Thus, we 676 change the key to standard key 677 For example: 678 'scf.XcType' -> 'scf_xctype' 679 """ 680 if isinstance(key, str): 681 return key.lower().replace('.', '_') 682 elif isinstance(key, list): 683 return [k.lower().replace('.', '_') for k in key] 684 else: 685 return [k.lower().replace('.', '_') for k in key] 686 687 688def get_standard_parameters(parameters): 689 """ 690 Translate the OpenMX parameters to standard ASE parameters. For example, 691 692 scf.XcType -> xc 693 scf.maxIter -> maxiter 694 scf.energycutoff -> energy_cutoff 695 scf.Kgrid -> kpts 696 scf.EigenvalueSolver -> eigensolver 697 scf.SpinPolarization -> spinpol 698 scf.criterion -> convergence 699 scf.Electric.Field -> external 700 scf.Mixing.Type -> mixer 701 scf.system.charge -> charge 702 703 We followed GPAW schem. 704 """ 705 from ase.calculators.openmx import parameters as param 706 from ase.units import Bohr, Ha, Ry, fs, m, s 707 units = param.unit_dat_keywords 708 standard_parameters = {} 709 standard_units = {'eV': 1, 'Ha': Ha, 'Ry': Ry, 'Bohr': Bohr, 'fs': fs, 710 'K': 1, 'GV / m': 1e9/1.6e-19 / m, 'Ha/Bohr': Ha/Bohr, 711 'm/s': m/s, '_amu': 1, 'Tesla': 1} 712 translated_parameters = { 713 'scf.XcType': 'xc', 714 'scf.maxIter': 'maxiter', 715 'scf.energycutoff': 'energy_cutoff', 716 'scf.Kgrid': 'kpts', 717 'scf.EigenvalueSolver': 'eigensolver', 718 'scf.SpinPolarization': 'spinpol', 719 'scf.criterion': 'convergence', 720 'scf.Electric.Field': 'external', 721 'scf.Mixing.Type': 'mixer', 722 'scf.system.charge': 'charge' 723 } 724 725 for key in parameters.keys(): 726 for openmx_key in translated_parameters.keys(): 727 if key == get_standard_key(openmx_key): 728 standard_key = translated_parameters[openmx_key] 729 unit = standard_units.get(units.get(openmx_key), 1) 730 standard_parameters[standard_key] = parameters[key] * unit 731 standard_parameters['spinpol'] = parameters.get('scf_spinpolarization') 732 return standard_parameters 733 734 735def get_atomic_formula(out_data=None, log_data=None, restart_data=None, 736 scfout_data=None, dat_data=None): 737 """_formula'. 738 OpenMX results gives following information. Since, we should pick one 739 between position/scaled_position, scaled_positions are suppressed by 740 default. We use input value of position. Not the position after 741 calculation. It is temporal. 742 743 Atoms.SpeciesAndCoordinate -> symbols 744 Atoms.SpeciesAndCoordinate -> positions 745 Atoms.UnitVectors -> cell 746 scaled_positions -> scaled_positions 747 If `positions` and `scaled_positions` are both given, this key deleted 748 magmoms -> magmoms, Single value for each atom or three numbers for each 749 atom for non-collinear calculations. 750 """ 751 atomic_formula = {} 752 parameters = {'symbols': list, 'positions': list, 'scaled_positions': list, 753 'magmoms': list, 'cell': list} 754 datas = [out_data, log_data, restart_data, scfout_data, dat_data] 755 atoms_unitvectors = None 756 atoms_spncrd_unit = 'ang' 757 atoms_unitvectors_unit = 'ang' 758 for data in datas: 759 # positions unit save 760 if 'atoms_speciesandcoordinates_unit' in data: 761 atoms_spncrd_unit = data['atoms_speciesandcoordinates_unit'] 762 # cell unit save 763 if 'atoms_unitvectors_unit' in data: 764 atoms_unitvectors_unit = data['atoms_unitvectors_unit'] 765 # symbols, positions or scaled_positions 766 if 'atoms_speciesandcoordinates' in data: 767 atoms_spncrd = data['atoms_speciesandcoordinates'] 768 # cell 769 if 'atoms_unitvectors' in data: 770 atoms_unitvectors = data['atoms_unitvectors'] 771 # pbc 772 if 'scf_eigenvaluesolver' in data: 773 scf_eigenvaluesolver = data['scf_eigenvaluesolver'] 774 # ??? 775 for openmx_keyword in data.keys(): 776 for standard_keyword in parameters.keys(): 777 if openmx_keyword == standard_keyword: 778 atomic_formula[standard_keyword] = data[openmx_keyword] 779 780 atomic_formula['symbols'] = [i[1] for i in atoms_spncrd] 781 782 openmx_spncrd_keyword = [[i[2], i[3], i[4]] for i in atoms_spncrd] 783 # Positions 784 positions_unit = atoms_spncrd_unit.lower() 785 positions = np.array(openmx_spncrd_keyword, dtype=float) 786 if positions_unit == 'ang': 787 atomic_formula['positions'] = positions 788 elif positions_unit == 'frac': 789 scaled_positions = np.array(openmx_spncrd_keyword, dtype=float) 790 atomic_formula['scaled_positions'] = scaled_positions 791 elif positions_unit == 'au': 792 positions = np.array(openmx_spncrd_keyword, dtype=float) * Bohr 793 atomic_formula['positions'] = positions 794 795 # If Cluster, pbc is False, else it is True 796 atomic_formula['pbc'] = scf_eigenvaluesolver.lower() != 'cluster' 797 798 # Cell Handling 799 if atoms_unitvectors is not None: 800 openmx_cell_keyword = atoms_unitvectors 801 cell = np.array(openmx_cell_keyword, dtype=float) 802 if atoms_unitvectors_unit.lower() == 'ang': 803 atomic_formula['cell'] = openmx_cell_keyword 804 elif atoms_unitvectors_unit.lower() == 'au': 805 atomic_formula['cell'] = cell * Bohr 806 807 # If `positions` and `scaled_positions` are both given, delete `scaled_..` 808 if atomic_formula.get('scaled_positions') is not None and \ 809 atomic_formula.get('positions') is not None: 810 del atomic_formula['scaled_positions'] 811 return atomic_formula 812 813 814def get_results(out_data=None, log_data=None, restart_data=None, 815 scfout_data=None, dat_data=None, band_data=None): 816 """ 817 From the gien data sets, construct the dictionary 'results' and return it' 818 OpenMX version 3.8 can yield following properties 819 free_energy, Ha # Same value with energy 820 energy, Ha 821 energies, Ha 822 forces, Ha/Bohr 823 stress(after 3.8 only) Ha/Bohr**3 824 dipole Debye 825 read_chemical_potential Ha 826 magmoms muB ?? set to 1 827 magmom muB ?? set to 1 828 """ 829 from numpy import array as arr 830 results = {} 831 implemented_properties = {'free_energy': Ha, 'energy': Ha, 'energies': Ha, 832 'forces': Ha/Bohr, 'stress': Ha/Bohr**3, 833 'dipole': Debye, 'chemical_potential': Ha, 834 'magmom': 1, 'magmoms': 1, 'eigenvalues': Ha} 835 data = [out_data, log_data, restart_data, scfout_data, dat_data, band_data] 836 for datum in data: 837 for key in datum.keys(): 838 for property in implemented_properties.keys(): 839 if key == property: 840 results[key] = arr(datum[key])*implemented_properties[key] 841 return results 842 843 844def get_file_name(extension='.out', filename=None, absolute_directory=True): 845 directory, prefix = os.path.split(filename) 846 if directory == '': 847 directory = os.curdir 848 if absolute_directory: 849 return os.path.abspath(directory + '/' + prefix + extension) 850 else: 851 return os.path.basename(directory + '/' + prefix + extension) 852