# flake8: noqa """ The ASE Calculator for OpenMX : Python interface to the software package for nano-scale material simulations based on density functional theories. Copyright (C) 2018 JaeHwan Shim and JaeJun Yu This program is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 2.1 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with ASE. If not, see . """ # from ase.calculators import SinglePointDFTCalculator import os import struct import numpy as np from ase.units import Ha, Bohr, Debye def read_openmx(filename=None, debug=False): from ase.calculators.openmx import OpenMX from ase import Atoms """ Read results from typical OpenMX output files and returns the atom object In default mode, it reads every implementd properties we could get from the files. Unlike previous version, we read the information based on file. previous results will be eraised unless previous results are written in the next calculation results. Read the 'LABEL.log' file seems redundant. Because all the information should already be written in '.out' file. However, in the version 3.8.3, stress tensor are not written in the '.out' file. It only contained in the '.log' file. So... I implented reading '.log' file method """ log_data = read_file(get_file_name('.log', filename), debug=debug) restart_data = read_file(get_file_name('.dat#', filename), debug=debug) dat_data = read_file(get_file_name('.dat', filename), debug=debug) out_data = read_file(get_file_name('.out', filename), debug=debug) scfout_data = read_scfout_file(get_file_name('.scfout', filename)) band_data = read_band_file(get_file_name('.Band', filename)) # dos_data = read_dos_file(get_file_name('.Dos.val', filename)) """ First, we get every data we could get from the all results files. And then, reform the data to fit to data structure of Atom object. While doing this, Fix the unit to ASE format. """ parameters = get_parameters(out_data=out_data, log_data=log_data, restart_data=restart_data, dat_data=dat_data, scfout_data=scfout_data, band_data=band_data) atomic_formula = get_atomic_formula(out_data=out_data, log_data=log_data, restart_data=restart_data, scfout_data=scfout_data, dat_data=dat_data) results = get_results(out_data=out_data, log_data=log_data, restart_data=restart_data, scfout_data=scfout_data, dat_data=dat_data, band_data=band_data) atoms = Atoms(**atomic_formula) atoms.calc = OpenMX(**parameters) atoms.calc.results = results return atoms def read_file(filename, debug=False): """ Read the 'LABEL.out' file. Using 'parameters.py', we read every 'allowed_ dat' dictionory. while reading a file, if one find the key matcheds That 'patters', which indicates the property we want is written, it will returns the pair value of that key. For example, example will be written later """ from ase.calculators.openmx import parameters as param if not os.path.isfile(filename): return {} param_keys = ['integer_keys', 'float_keys', 'string_keys', 'bool_keys', 'list_int_keys', 'list_float_keys', 'list_bool_keys', 'tuple_integer_keys', 'tuple_float_keys', 'tuple_float_keys'] patterns = { 'Stress tensor': ('stress', read_stress_tensor), 'Dipole moment': ('dipole', read_dipole), 'Fractional coordinates of': ('scaled_positions', read_scaled_positions), 'Utot.': ('energy', read_energy), 'energies in': ('energies', read_energies), 'Chemical Potential': ('chemical_potential', read_chemical_potential), '>> str(rn(line, 1)) LDA >>> line = f.readline() >>> int(rn(line, 3)) 4 """ return line.split()[-n] def read_tuple_integer(line): return tuple([int(x) for x in line.split()[-3:]]) def read_tuple_float(line): return tuple([float(x) for x in line.split()[-3:]]) def read_integer(line): return int(rn(line)) def read_float(line): return float(rn(line)) def read_string(line): return str(rn(line)) def read_bool(line): bool = str(rn(line)).lower() if bool == 'on': return True elif bool == 'off': return False else: return None def read_list_int(line): return [int(x) for x in line.split()[1:]] def read_list_float(line): return [float(x) for x in line.split()[1:]] def read_list_bool(line): return [read_bool(x) for x in line.split()[1:]] def read_matrix(line, key, f): matrix = [] line = f.readline() while key not in line: matrix.append(line.split()) line = f.readline() return matrix def read_stress_tensor(line, f, debug=None): f.readline() # passing empty line f.readline() line = f.readline() xx, xy, xz = read_tuple_float(line) line = f.readline() yx, yy, yz = read_tuple_float(line) line = f.readline() zx, zy, zz = read_tuple_float(line) stress = [xx, yy, zz, (zy + yz)/2, (zx + xz)/2, (yx + xy)/2] return stress def read_magmoms_and_total_magmom(line, f, debug=None): total_magmom = read_float(line) f.readline() # Skip empty lines f.readline() line = f.readline() magmoms = [] while not(line == '' or line.isspace()): magmoms.append(read_float(line)) line = f.readline() return magmoms, total_magmom def read_energy(line, f, debug=None): # It has Hartree unit yet return read_float(line) def read_energies(line, f, debug=None): line = f.readline() if '***' in line: point = 7 # Version 3.8 else: point = 16 # Version 3.9 for i in range(point): f.readline() line = f.readline() energies = [] while not(line == '' or line.isspace()): energies.append(float(line.split()[2])) line = f.readline() return energies def read_eigenvalues(line, f, debug=False): """ Read the Eigenvalues in the `.out` file and returns the eigenvalue First, it assumes system have two spins and start reading until it reaches the end('*****...'). eigenvalues[spin][kpoint][nbands] For symmetry reason, `.out` file prints the eigenvalues at the half of the K points. Thus, we have to fill up the rest of the half. However, if the calculation was conducted only on the gamma point, it will raise the 'gamma_flag' as true and it will returns the original samples. """ def prind(line): if debug: print(line) if 'Hartree' in line: return None prind("Read eigenvalue output") current_line = f.tell() f.seek(0) # Seek for the kgrid information while line != '': line = f.readline().lower() if 'scf.kgrid' in line: break f.seek(current_line) # Retrun to the original position kgrid = read_tuple_integer(line) line = f.readline() line = f.readline() if kgrid != (): prind('Non-Gamma point calculation') prind('scf.Kgrid is %d, %d, %d' % kgrid) gamma_flag = False f.seek(f.tell()+57) else: prind('Gamma point calculation') gamma_flag = True eigenvalues = [] eigenvalues.append([]) eigenvalues.append([]) # Assume two spins i = 0 while 'Mulliken' not in line: line = f.readline() prind(line) eigenvalues[0].append([]) eigenvalues[1].append([]) while not (line == '' or line.isspace()): eigenvalues[0][i].append(float(rn(line, 2))) eigenvalues[1][i].append(float(rn(line, 1))) line = f.readline() prind(line) i += 1 f.readline() f.readline() line = f.readline() prind(line) if gamma_flag: return np.asarray(eigenvalues) eigen_half = np.asarray(eigenvalues) prind(eigen_half) # Fill up the half spin, half_kpts, bands = eigen_half.shape even_odd = np.array(kgrid).prod() % 2 eigen_values = np.zeros((spin, half_kpts*2-even_odd, bands)) for i in range(half_kpts): eigen_values[0, i] = eigen_half[0, i, :] eigen_values[1, i] = eigen_half[1, i, :] eigen_values[0, 2*half_kpts-1-i-even_odd] = eigen_half[0, i, :] eigen_values[1, 2*half_kpts-1-i-even_odd] = eigen_half[1, i, :] return eigen_values def read_forces(line, f, debug=None): # It has Hartree per Bohr unit yet forces = [] f.readline() # Skip Empty line line = f.readline() while 'coordinates.forces>' not in line: forces.append(read_tuple_float(line)) line = f.readline() return np.array(forces) def read_dipole(line, f, debug=None): dipole = [] while 'Total' not in line: line = f.readline() dipole.append(read_tuple_float(line)) return dipole def read_scaled_positions(line, f, debug=None): scaled_positions = [] f.readline() # Skip Empty lines f.readline() f.readline() line = f.readline() while not(line == '' or line.isspace()): # Detect empty line scaled_positions.append(read_tuple_float(line)) line = f.readline() return scaled_positions def read_chemical_potential(line, f, debug=None): return read_float(line) def get_parameters(out_data=None, log_data=None, restart_data=None, scfout_data=None, dat_data=None, band_data=None): """ From the given data sets, construct the dictionary 'parameters'. If data is in the paramerters, it will save it. """ from ase.calculators.openmx import parameters as param scaned_data = [dat_data, out_data, log_data, restart_data, scfout_data, band_data] openmx_keywords = [param.tuple_integer_keys, param.tuple_float_keys, param.tuple_bool_keys, param.integer_keys, param.float_keys, param.string_keys, param.bool_keys, param.list_int_keys, param.list_bool_keys, param.list_float_keys, param.matrix_keys] parameters = {} for scaned_datum in scaned_data: for scaned_key in scaned_datum.keys(): for openmx_keyword in openmx_keywords: if scaned_key in get_standard_key(openmx_keyword): parameters[scaned_key] = scaned_datum[scaned_key] continue translated_parameters = get_standard_parameters(parameters) parameters.update(translated_parameters) return {k: v for k, v in parameters.items() if v is not None} def get_standard_key(key): """ Standard ASE parameter format is to USE unerbar(_) instead of dot(.). Also, It is recommended to use lower case alphabet letter. Not Upper. Thus, we change the key to standard key For example: 'scf.XcType' -> 'scf_xctype' """ if isinstance(key, str): return key.lower().replace('.', '_') elif isinstance(key, list): return [k.lower().replace('.', '_') for k in key] else: return [k.lower().replace('.', '_') for k in key] def get_standard_parameters(parameters): """ Translate the OpenMX parameters to standard ASE parameters. For example, scf.XcType -> xc scf.maxIter -> maxiter scf.energycutoff -> energy_cutoff scf.Kgrid -> kpts scf.EigenvalueSolver -> eigensolver scf.SpinPolarization -> spinpol scf.criterion -> convergence scf.Electric.Field -> external scf.Mixing.Type -> mixer scf.system.charge -> charge We followed GPAW schem. """ from ase.calculators.openmx import parameters as param from ase.units import Bohr, Ha, Ry, fs, m, s units = param.unit_dat_keywords standard_parameters = {} standard_units = {'eV': 1, 'Ha': Ha, 'Ry': Ry, 'Bohr': Bohr, 'fs': fs, 'K': 1, 'GV / m': 1e9/1.6e-19 / m, 'Ha/Bohr': Ha/Bohr, 'm/s': m/s, '_amu': 1, 'Tesla': 1} translated_parameters = { 'scf.XcType': 'xc', 'scf.maxIter': 'maxiter', 'scf.energycutoff': 'energy_cutoff', 'scf.Kgrid': 'kpts', 'scf.EigenvalueSolver': 'eigensolver', 'scf.SpinPolarization': 'spinpol', 'scf.criterion': 'convergence', 'scf.Electric.Field': 'external', 'scf.Mixing.Type': 'mixer', 'scf.system.charge': 'charge' } for key in parameters.keys(): for openmx_key in translated_parameters.keys(): if key == get_standard_key(openmx_key): standard_key = translated_parameters[openmx_key] unit = standard_units.get(units.get(openmx_key), 1) standard_parameters[standard_key] = parameters[key] * unit standard_parameters['spinpol'] = parameters.get('scf_spinpolarization') return standard_parameters def get_atomic_formula(out_data=None, log_data=None, restart_data=None, scfout_data=None, dat_data=None): """_formula'. OpenMX results gives following information. Since, we should pick one between position/scaled_position, scaled_positions are suppressed by default. We use input value of position. Not the position after calculation. It is temporal. Atoms.SpeciesAndCoordinate -> symbols Atoms.SpeciesAndCoordinate -> positions Atoms.UnitVectors -> cell scaled_positions -> scaled_positions If `positions` and `scaled_positions` are both given, this key deleted magmoms -> magmoms, Single value for each atom or three numbers for each atom for non-collinear calculations. """ atomic_formula = {} parameters = {'symbols': list, 'positions': list, 'scaled_positions': list, 'magmoms': list, 'cell': list} datas = [out_data, log_data, restart_data, scfout_data, dat_data] atoms_unitvectors = None atoms_spncrd_unit = 'ang' atoms_unitvectors_unit = 'ang' for data in datas: # positions unit save if 'atoms_speciesandcoordinates_unit' in data: atoms_spncrd_unit = data['atoms_speciesandcoordinates_unit'] # cell unit save if 'atoms_unitvectors_unit' in data: atoms_unitvectors_unit = data['atoms_unitvectors_unit'] # symbols, positions or scaled_positions if 'atoms_speciesandcoordinates' in data: atoms_spncrd = data['atoms_speciesandcoordinates'] # cell if 'atoms_unitvectors' in data: atoms_unitvectors = data['atoms_unitvectors'] # pbc if 'scf_eigenvaluesolver' in data: scf_eigenvaluesolver = data['scf_eigenvaluesolver'] # ??? for openmx_keyword in data.keys(): for standard_keyword in parameters.keys(): if openmx_keyword == standard_keyword: atomic_formula[standard_keyword] = data[openmx_keyword] atomic_formula['symbols'] = [i[1] for i in atoms_spncrd] openmx_spncrd_keyword = [[i[2], i[3], i[4]] for i in atoms_spncrd] # Positions positions_unit = atoms_spncrd_unit.lower() positions = np.array(openmx_spncrd_keyword, dtype=float) if positions_unit == 'ang': atomic_formula['positions'] = positions elif positions_unit == 'frac': scaled_positions = np.array(openmx_spncrd_keyword, dtype=float) atomic_formula['scaled_positions'] = scaled_positions elif positions_unit == 'au': positions = np.array(openmx_spncrd_keyword, dtype=float) * Bohr atomic_formula['positions'] = positions # If Cluster, pbc is False, else it is True atomic_formula['pbc'] = scf_eigenvaluesolver.lower() != 'cluster' # Cell Handling if atoms_unitvectors is not None: openmx_cell_keyword = atoms_unitvectors cell = np.array(openmx_cell_keyword, dtype=float) if atoms_unitvectors_unit.lower() == 'ang': atomic_formula['cell'] = openmx_cell_keyword elif atoms_unitvectors_unit.lower() == 'au': atomic_formula['cell'] = cell * Bohr # If `positions` and `scaled_positions` are both given, delete `scaled_..` if atomic_formula.get('scaled_positions') is not None and \ atomic_formula.get('positions') is not None: del atomic_formula['scaled_positions'] return atomic_formula def get_results(out_data=None, log_data=None, restart_data=None, scfout_data=None, dat_data=None, band_data=None): """ From the gien data sets, construct the dictionary 'results' and return it' OpenMX version 3.8 can yield following properties free_energy, Ha # Same value with energy energy, Ha energies, Ha forces, Ha/Bohr stress(after 3.8 only) Ha/Bohr**3 dipole Debye read_chemical_potential Ha magmoms muB ?? set to 1 magmom muB ?? set to 1 """ from numpy import array as arr results = {} implemented_properties = {'free_energy': Ha, 'energy': Ha, 'energies': Ha, 'forces': Ha/Bohr, 'stress': Ha/Bohr**3, 'dipole': Debye, 'chemical_potential': Ha, 'magmom': 1, 'magmoms': 1, 'eigenvalues': Ha} data = [out_data, log_data, restart_data, scfout_data, dat_data, band_data] for datum in data: for key in datum.keys(): for property in implemented_properties.keys(): if key == property: results[key] = arr(datum[key])*implemented_properties[key] return results def get_file_name(extension='.out', filename=None, absolute_directory=True): directory, prefix = os.path.split(filename) if directory == '': directory = os.curdir if absolute_directory: return os.path.abspath(directory + '/' + prefix + extension) else: return os.path.basename(directory + '/' + prefix + extension)