From 6b69257f30a2ee6ece8929ba8d12dc6cfa654425 Mon Sep 17 00:00:00 2001 From: Alvin Noe Ladines <ladinesalvinnoe@gmail.com> Date: Fri, 22 Jan 2021 03:28:06 +0100 Subject: [PATCH] Fix reusing of parser (issue #485) --- elasticparser/elastic_parser.py | 544 ++++++++++++++++++++++++++-- elasticparser/elastic_properties.py | 466 ------------------------ 2 files changed, 507 insertions(+), 503 deletions(-) delete mode 100644 elasticparser/elastic_properties.py diff --git a/elasticparser/elastic_parser.py b/elasticparser/elastic_parser.py index ceb78d4..4f590f1 100644 --- a/elasticparser/elastic_parser.py +++ b/elasticparser/elastic_parser.py @@ -1,39 +1,479 @@ import os import numpy as np import logging +import pint +from ase import Atoms +from nomad.parsing.parser import FairdiParser +from nomad.parsing.file_parser import Quantity, TextParser from nomad.datamodel.metainfo.common_dft import Run, System, SingleConfigurationCalculation,\ Method, CalculationToCalculationRefs, Workflow, Elastic -from nomad.parsing.parser import FairdiParser from elasticparser.metainfo.elastic import x_elastic_section_strain_diagrams,\ x_elastic_section_fitting_parameters -from elasticparser.elastic_properties import ElasticProperties from .metainfo import m_env -class ElasticParser(FairdiParser): +class InfoParser(TextParser): def __init__(self): - super().__init__( - name='parsers/elastic', code_name='elastic', code_homepage='http://exciting-code.org/elastic', - mainfile_contents_re=r'\s*Order of elastic constants\s*=\s*[0-9]+\s*', - mainfile_name_re=(r'.*/INFO_ElaStic')) + super().__init__(None) + + def init_quantities(self): + self._quantities = [ + Quantity( + 'order', r'\s*Order of elastic constants\s*=\s*([0-9]+)', repeats=False, + dtype=int), + Quantity( + 'calculation_method', r'\s*Method of calculation\s*=\s*([-a-zA-Z]+)\s*', + repeats=False), + Quantity( + 'code_name', r'\s*DFT code name\s*=\s*([-a-zA-Z]+)', repeats=False), + Quantity( + 'space_group_number', r'\s*Space-group number\s*=\s*([0-9]+)', repeats=False), + Quantity( + 'equilibrium_volume', r'\s*Volume of equilibrium unit cell\s*=\s*([0-9.]+)\s*', + unit='angstrom ** 3'), + Quantity( + 'max_strain', r'\s*Maximum Lagrangian strain\s*=\s*([0-9.]+)', repeats=False), + Quantity( + 'n_strains', r'\s*Number of distorted structures\s*=\s*([0-9]+)', repeats=False)] + + +class StructureParser(TextParser): + def __init__(self): + super().__init__(None) + + def init_quantities(self): + def get_sym_pos(val): + val = val.strip().replace('\n', '').split() + sym = [] + pos = [] + for i in range(0, len(val), 4): + sym.append(val[i + 3].strip()) + pos.append([float(val[j]) for j in range(i, i + 3)]) + sym_pos = dict(symbols=sym, positions=pos) + return sym_pos + + self._quantities = [ + Quantity( + 'cellpar', r'a\s*b\s*c\n([\d\.\s]+)\n\s*alpha\s*beta\s*gamma\n([\d\.\s]+)\n+', + repeats=False), + Quantity( + 'sym_pos', r'Atom positions:\n\n([\s\d\.A-Za-z]+)\n\n', + str_operation=get_sym_pos, repeats=False, convert=False)] + + +class DistortedParametersParser(TextParser): + def __init__(self): + super().__init__(None) - self._metainfo_env = m_env - self.properties = ElasticProperties() + def init_quantities(self): + self._quantities = [Quantity( + 'deformation', r'Lagrangian strain\s*=\s*\(([eta\s\d\.,]+)\)', + str_operation=lambda x: x.replace(',', '').split(), repeats=True, dtype=str)] + + +class FitParser(TextParser): + def __init__(self): + super().__init__(None) + + def init_quantities(self): + + def split_eta_val(val): + order, val = val.strip().split(' order fit.') + val = [float(v) for v in val.strip().split()] + return order, val[0::2], val[1::2] + + self._quantities = [Quantity( + 'fit', r'(\w+ order fit\.\n[\d.\s\neE\-\+]+)\n', repeats=True, convert=False, + str_operation=split_eta_val)] + + +class ElasticConstant2Parser(TextParser): + def __init__(self): + super().__init__(None) + + def init_quantities(self): + self._quantities = [ + Quantity( + 'voigt', r'Symmetry[\s\S]+\n\s*\n([C\d\s\n\(\)\-\+\/\*]+)\n', + shape=(6, 6), dtype=str, repeats=False), + Quantity( + 'elastic_constant', r'Elastic constant[\s\S]+in GPa\s*:\s*\n\n([\-\d\.\s\n]+)\n', + shape=(6, 6), dtype=float, unit='GPa', repeats=False), + Quantity( + 'compliance', r'Elastic compliance[\s\S]+in 1/GPa\s*:\s*\n\n([\-\d\.\s\n]+)\n', + shape=(6, 6), dtype=float, unit='1/GPa', repeats=False)] + + def str_to_modulus(val_in): + val_in = val_in.strip().split() + key = val_in[0] + unit = val_in[-1] if len(val_in) == 3 else None + val = float(val_in[1]) + val = pint.Quantity(val, unit) if unit is not None else val + return key, val + + self._quantities.append(Quantity( + 'modulus', r',\s*(\w+)\s*=\s*([\-\+\w\. ]+?)\n', str_operation=str_to_modulus, + repeats=True)) + + self._quantities.append(Quantity( + 'eigenvalues', + r'Eigenvalues of elastic constant \(stiffness\) matrix:\s*\n+([\-\d\.\n\s]+)\n', + unit='GPa', repeats=False)) + + +class ElasticConstant3Parser(TextParser): + def __init__(self): + super().__init__(None) + + def init_quantities(self): + def arrange_matrix(val): + val = val.strip().split('\n') + matrix = [v.strip().split() for v in val if v.strip()] + matrix = np.array(matrix).reshape((12, 18)) + arranged = [] + for i in range(2): + for j in range(3): + arranged.append( + matrix[i * 6: (i + 1) * 6, j * 6: (j + 1) * 6].tolist()) + return arranged + + self._quantities = [ + Quantity( + 'elastic_constant', r'\%\s*\n([\s0-6A-L]*)[\n\s\%1-6\-ij]*([\s0-6A-L]*)\n', + str_operation=arrange_matrix, dtype=str, repeats=False, convert=False), + Quantity( + 'cijk', r'(C\d\d\d)\s*=\s*([\-\d\.]+)\s*GPa', repeats=True, convert=False)] + + +class ElasticParserInterface: + def __init__(self): + self._mainfile = None + self.logger = None + self._deform_dirs = None + self._deform_dir_prefix = 'Dst' + self.info = InfoParser() + self.structure = StructureParser() + self.distorted_parameters = DistortedParametersParser() + self.fit = FitParser() + self.elastic_constant_2 = ElasticConstant2Parser() + self.elastic_constant_3 = ElasticConstant3Parser() + + @property + def deformation_dirs(self): + if self._deform_dirs is None: + dirs = os.listdir(self.maindir) + self._deform_dirs = [ + os.path.join(self.maindir, d) for d in dirs if d.startswith(self._deform_dir_prefix)] + + return self._deform_dirs + + def get_references_to_calculations(self): + def output_file(dirname): + code = self.info.get('code_name', '').lower() + if code == 'exciting': + return os.path.join(dirname, 'INFO.OUT') + elif code == 'wien': + return os.path.join(dirname, '%s_Converged.scf' % os.path.basename(dirname)) + elif code == 'quantum': + return os.path.join(dirname, '%s.out' % os.path.basename(dirname)) + else: + return None + + references = [] + for deform_dir in self.deformation_dirs: + sub_dirs = os.listdir(deform_dir) + for sub_dir in sub_dirs: + calc_dir = os.path.join(deform_dir, sub_dir) + out_file = output_file(calc_dir) + if out_file is not None and os.path.isfile(out_file): + references.append(out_file) + + return references + + def get_structure_info(self): + path = os.path.join(self.maindir, 'sgroup.out') + if not os.path.isfile(path): + return + + self.structure.mainfile = path + + cellpar = self.structure.get('cellpar', None) + sym_pos = self.structure.get('sym_pos', {}) + + sym = sym_pos.get('symbols', None) + pos = sym_pos.get('positions', None) + + if cellpar is None or sym is None or pos is None: + return + + structure = Atoms(cell=cellpar, scaled_positions=pos, symbols=sym, pbc=True) + + positions = structure.get_positions() + positions = pint.Quantity(positions, 'angstrom') + cell = structure.get_cell() + cell = pint.Quantity(cell, 'angstrom') + + return sym, positions, cell + + def get_strain_energy(self): + strains, energies = [], [] + + for deform_dir in self.deformation_dirs: + filenames = [d for d in os.listdir(deform_dir) if d.endswith('Energy.dat')] + if not filenames: + continue + + path = os.path.join(deform_dir, filenames[-1]) + data = np.loadtxt(path).T + strains.append(list(data[0])) + # the peculiarity of the x_elastic_strain_diagram_values metainfo that it does + # not have the energy unit + energies.append(list(pint.Quantity(data[1], 'hartree').to('J').magnitude)) + + return strains, energies + + def get_strain_stress(self): + strains = {'Lagrangian-stress': [], 'Physical-stress': []} + stresses = {'Lagrangian-stress': [], 'Physical-stress': []} + + for deform_dir in self.deformation_dirs: + filenames = [d for d in os.listdir(deform_dir) if d.endswith('stress.dat')] + + for filename in filenames: + path = os.path.join(deform_dir, filename) + if not os.path.isfile(path): + continue + + with open(path) as f: + lines = f.readlines() + + strain, stress = [], [] + for line in lines: + val = line.strip().split() + if not val[0].strip().replace('.', '').isdecimal(): + continue + + strain.append(float(val[0])) + stress.append([float(v) for v in val[1:7]]) + + stype = filename.rstrip('.dat').split('_')[-1] + strains[stype].append(strain) + stresses[stype].append(stress) + + return strains, stresses + + def get_deformation_types(self): + path = os.path.join(self.maindir, 'Distorted_Parameters') + self.distorted_parameters.mainfile = path + return self.distorted_parameters.get('deformation') + + def _get_fit(self, path_dir, file_ext): + path_dir = os.path.join(self.maindir, path_dir) + + if not os.path.isdir(path_dir): + return + + paths = [p for p in os.listdir(path_dir) if p.endswith(file_ext)] + paths.sort() + + if not paths: + return + + eta, val = {}, {} + for path in paths: + self.fit.mainfile = os.path.join(path_dir, path) + fit_results = self.fit.get('fit', []) + for result in fit_results: + eta.setdefault(result[0], []) + val.setdefault(result[0], []) + eta[result[0]].append(result[1]) + val[result[0]].append(result[2]) + + return eta, val + + def get_energy_fit(self): + energy_fit = dict() + + for file_ext in ['d2E.dat', 'd3E.dat', 'ddE.dat']: + result = self._get_fit('Energy-vs-Strain', file_ext) + if result is None: + continue + + result = list(result) + result[1] = { + key: pint.Quantity( + val, 'GPa').to('Pa').magnitude for key, val in result[1].items()} + energy_fit['d2e'] = result + + result = self._get_fit('Energy-vs-Strain', 'CVe.dat') + if result is not None: + result = list(result) + result[1] = { + key: pint.Quantity( + val, 'hartree').to('J').magnitude for key, val in result[1].items()} + energy_fit['cross-validation'] = result + + return energy_fit + + def get_stress_fit(self): + stress_fit = dict() + stress_fit['dtn'] = [[]] * 6 + stress_fit['cross-validation'] = [[]] * 6 + + for strain_index in range(1, 7): + result = self._get_fit('Stress-vs-Strain', '%d_dS.dat' % strain_index) + if result is not None: + result[1] = {key: pint.Quantity(val, 'GPa') for key, val in result[1].items()} + stress_fit['dtn'][strain_index - 1] = result + + result = self._get_fit('Stress-vs-Strain', '%d_CVe.dat' % strain_index) + if result is not None: + result[1] = {key: pint.Quantity(val, 'hartree') for key, val in result[1].items()} + stress_fit['cross-validation'][strain_index - 1] = result + + return stress_fit + + def get_input(self): + paths = os.listdir(self.maindir) + path = None + order = self.info.get('order', 2) + for p in paths: + if 'ElaStic_' in p and p.endswith('.in') and str(order) in p: + path = p + break + + if path is None: + return + + calc_method = self.info.get('calculation_method') + + eta_ec = [] + fit_ec = [] + + def _is_number(var): + try: + float(var) + return True + except Exception: + return False + + path = os.path.join(self.maindir, path) + + with open(path) as f: + while True: + line = f.readline() + if not line: + break + + if calc_method.lower() == 'energy': + _, eta, fit = line.strip().split() + eta_ec.append(float(eta)) + fit_ec.append(int(fit)) + + elif calc_method.lower() == 'stress': + val = line.strip().split() + if not _is_number(val[0]): + eta_ec.append([float(val[i + 1]) for i in range(6)]) + else: + fit_ec.append([int(val[i]) for i in range(6)]) + + else: + pass + + return eta_ec, fit_ec + + def get_elastic_constants_order2(self): + path = os.path.join(self.maindir, 'ElaStic_2nd.out') + + self.elastic_constant_2.mainfile = path + + matrices = dict() + for key in ['voigt', 'elastic_constant', 'compliance']: + val = self.elastic_constant_2.get(key, None) + if val is not None: + matrices[key] = val + + moduli = dict() + for modulus in self.elastic_constant_2.get('modulus', []): + moduli[modulus[0]] = modulus[1] + + eigenvalues = self.elastic_constant_2.get('eigenvalues') + + return matrices, moduli, eigenvalues + + def get_elastic_constants_order3(self): + path = os.path.join(self.maindir, 'ElaStic_3rd.out') + self.elastic_constant_3.mainfile = path + + elastic_constant_str = self.elastic_constant_3.get('elastic_constant') + if elastic_constant_str is None: + return + + cijk = dict() + for element in self.elastic_constant_3.get('cijk', []): + cijk[element[0]] = float(element[1]) + + # formulas for the coefficients + coeff_A = cijk.get('C111', 0) + cijk.get('C112', 0) - cijk.get('C222', 0) + coeff_B = -(cijk.get('C115', 0) + 3 * cijk.get('C125', 0)) / 2 + coeff_C = (cijk.get('C114', 0) + 3 * cijk.get('C124', 0)) / 2 + coeff_D = -(2 * cijk.get('C111', 0) + cijk.get('C112', 0) - 3 * cijk.get('C222', 0)) / 4 + coeff_E = -cijk.get('C114', 0) - 2 * cijk.get('C124', 0) + coeff_F = -cijk.get('C115', 0) - 2 * cijk.get('C125', 0) + coeff_G = -(cijk.get('C115', 0) - cijk.get('C125', 0)) / 2 + coeff_H = (cijk.get('C114', 0) - cijk.get('C124', 0)) / 2 + coeff_I = (2 * cijk.get('C111', 0) - cijk.get('C112', 0) - cijk.get('C222', 0)) / 4 + coeff_J = (cijk.get('C113', 0) - cijk.get('C123', 0)) / 2 + coeff_K = -(cijk.get('C144', 0) - cijk.get('C155', 0)) / 2 + + space_group_number = self.info.get('space_group_number') + + if space_group_number <= 148: # rhombohedral II + coefficients = dict( + A=coeff_A, B=coeff_B, C=coeff_C, D=coeff_D, E=coeff_E, F=coeff_F, G=coeff_G, + H=coeff_H, I=coeff_I, J=coeff_J, K=coeff_K) + elif space_group_number <= 167: # rhombohedral I + coefficients = dict( + A=coeff_A, B=coeff_C, C=coeff_D, D=coeff_E, E=coeff_H, F=coeff_I, G=coeff_J, + H=coeff_K) + elif space_group_number <= 176: # hexagonal II + coefficients = dict( + A=coeff_A, B=coeff_D, C=coeff_I, D=coeff_J, E=coeff_K) + elif space_group_number <= 194: # hexagonal I + coefficients = dict( + A=coeff_A, B=coeff_D, C=coeff_I, D=coeff_J, E=coeff_K) + else: + coefficients = dict() + + # assign values to the elastic constant matrix from the independent ec and coefficients + elastic_constant = np.zeros((6, 6, 6)) + for i in range(6): + for j in range(6): + for k in range(6): + val = elastic_constant_str[i][j][k] + if val == '0': + elastic_constant[i][j][k] = 0 + elif val.isdigit(): + elastic_constant[i][j][k] = cijk['C%s' % val] + else: + elastic_constant[i][j][k] = coefficients.get(val, 0) + + return elastic_constant def parse_strain(self): sec_scc = self.archive.section_run[-1].section_single_configuration_calculation[-1] - method = self.properties.info['calculation_method'].lower() + method = self.info['calculation_method'].lower() if method == 'energy': - strain, energy = self.properties.get_strain_energy() + strain, energy = self.get_strain_energy() if not strain: self.logger.warn('Error getting strain and energy data') return - n_strains = self.properties.info['n_strains'] + n_strains = self.info['n_strains'] sec_strain_diagram = sec_scc.m_create(x_elastic_section_strain_diagrams) sec_strain_diagram.x_elastic_strain_diagram_type = 'energy' @@ -47,7 +487,7 @@ class ElasticParser(FairdiParser): '2nd': poly_fit_2, '3rd': poly_fit_2 - 1, '4th': poly_fit_2 - 1, '5th': poly_fit_2 - 2, '6th': poly_fit_2 - 2, '7th': poly_fit_2 - 3} - energy_fit = self.properties.get_energy_fit() + energy_fit = self.get_energy_fit() if not energy_fit: self.logger.warn('Error getting energy fit data') return @@ -62,7 +502,7 @@ class ElasticParser(FairdiParser): sec_strain_diagram.x_elastic_strain_diagram_values = energy_fit[diagram_type][1][fit_order] elif method == 'stress': - strain, stress = self.properties.get_strain_stress() + strain, stress = self.get_strain_stress() for diagram_type in ['Lagrangian-stress', 'Physical-stress']: strain_i = strain[diagram_type] stress_i = np.transpose(np.array(stress[diagram_type]), axes=(2, 0, 1)) @@ -77,7 +517,7 @@ class ElasticParser(FairdiParser): sec_strain_diagram.x_elastic_strain_diagram_eta_values = strain_i sec_strain_diagram.x_elastic_strain_diagram_values = stress_i[si] - stress_fit = self.properties.get_stress_fit() + stress_fit = self.get_stress_fit() for diagram_type in ['cross-validation', 'dtn']: if stress_fit.get(diagram_type, None) is None: continue @@ -95,10 +535,10 @@ class ElasticParser(FairdiParser): def parse_elastic_constant(self): sec_scc = self.archive.section_run[-1].section_single_configuration_calculation[-1] - order = self.properties.info['order'] + order = self.info['order'] if order == 2: - matrices, moduli, eigenvalues = self.properties.get_elastic_constants_order2() + matrices, moduli, eigenvalues = self.get_elastic_constants_order2() sec_scc.x_elastic_2nd_order_constants_notation_matrix = matrices['voigt'] sec_scc.x_elastic_2nd_order_constants_matrix = matrices['elastic_constant'] sec_scc.x_elastic_2nd_order_constants_compliance_matrix = matrices['compliance'] @@ -122,20 +562,29 @@ class ElasticParser(FairdiParser): sec_scc.x_elastic_eigenvalues = eigenvalues elif order == 3: - elastic_constant = self.properties.get_elastic_constants_order3() + elastic_constant = self.get_elastic_constants_order3() sec_scc.x_elastic_3rd_order_constants_matrix = elastic_constant - def _init_parsers(self): - self.properties.mainfile = self.filepath - self.properties.logger = self.logger + def init_parser(self): + self._deform_dirs = None + self.maindir = os.path.dirname(self.filepath) + self.info.mainfile = self.filepath + + def reuse_parser(self, parser): + self.info.quantities = parser.quantities + self.structure.quantities = parser.structure.quantities + self.distorted_parameters.quantities = parser.distorted_parameters.quantities + self.fit.quantities = parser.fit.quantities + self.elastic_constant_2.quantities = parser.elastic_constant_2.quantities + self.elastic_constant_3.quantities = parser.elastic_constant_3.quantities def parse(self, filepath, archive, logger): self.filepath = os.path.abspath(filepath) self.archive = archive self.logger = logger if logger is not None else logging - self._init_parsers() + self.init_parser() sec_run = self.archive.m_create(Run) @@ -144,35 +593,35 @@ class ElasticParser(FairdiParser): sec_system = sec_run.m_create(System) - symbols, positions, cell = self.properties.get_structure_info() - volume = self.properties.info['equilibrium_volume'] + symbols, positions, cell = self.get_structure_info() + volume = self.info['equilibrium_volume'] sec_system.atom_labels = symbols sec_system.atom_positions = positions sec_system.simulation_cell = cell sec_system.configuration_periodic_dimensions = [True, True, True] - sec_system.x_elastic_space_group_number = self.properties.info['space_group_number'] + sec_system.x_elastic_space_group_number = self.info['space_group_number'] sec_system.x_elastic_unit_cell_volume = volume sec_method = sec_run.m_create(Method) - sec_method.x_elastic_elastic_constant_order = self.properties.info['order'] - sec_method.x_elastic_calculation_method = self.properties.info['calculation_method'] - sec_method.x_elastic_code = self.properties.info['code_name'] - sec_method.x_elastic_max_lagrangian_strain = self.properties.info['max_strain'] - sec_method.x_elastic_number_of_distorted_structures = self.properties.info['n_strains'] + sec_method.x_elastic_elastic_constant_order = self.info['order'] + sec_method.x_elastic_calculation_method = self.info['calculation_method'] + sec_method.x_elastic_code = self.info['code_name'] + sec_method.x_elastic_max_lagrangian_strain = self.info['max_strain'] + sec_method.x_elastic_number_of_distorted_structures = self.info['n_strains'] - deformation_types = self.properties.get_deformation_types() + deformation_types = self.get_deformation_types() sec_method.x_elastic_deformation_types = deformation_types - sec_method.x_elastic_number_of_deformations = len(self.properties.deformation_dirs) + sec_method.x_elastic_number_of_deformations = len(self.deformation_dirs) - references = self.properties.get_references_to_calculations() + references = self.get_references_to_calculations() sec_scc = sec_run.m_create(SingleConfigurationCalculation) for reference in references: sec_calc_ref = sec_scc.m_create(CalculationToCalculationRefs) sec_calc_ref.calculation_to_calculation_external_url = reference sec_calc_ref.calculation_to_calculation_kind = 'source_calculation' - fit_input = self.properties.get_input() + fit_input = self.get_input() sec_fit_par = sec_method.m_create(x_elastic_section_fitting_parameters) sec_fit_par.x_elastic_fitting_parameters_eta = fit_input[0] sec_fit_par.x_elastic_fitting_parameters_polynomial_order = fit_input[1] @@ -187,6 +636,27 @@ class ElasticParser(FairdiParser): sec_workflow = self.archive.m_create(Workflow) sec_workflow.workflow_type = 'elastic' sec_elastic = sec_workflow.m_create(Elastic) - sec_elastic.elastic_calculation_method = self.properties.info['calculation_method'].lower() - sec_elastic.elastic_constants_order = self.properties.info['order'] - sec_elastic.strain_maximum = self.properties.info['max_strain'] + sec_elastic.elastic_calculation_method = self.info['calculation_method'].lower() + sec_elastic.elastic_constants_order = self.info['order'] + sec_elastic.strain_maximum = self.info['max_strain'] + + +class ElasticParser(FairdiParser): + def __init__(self): + super().__init__( + name='parsers/elastic', code_name='elastic', code_homepage='http://exciting-code.org/elastic', + mainfile_contents_re=r'\s*Order of elastic constants\s*=\s*[0-9]+\s*', + mainfile_name_re=(r'.*/INFO_ElaStic')) + + self._metainfo_env = m_env + self.parser = None + + def parse(self, filepath, archive, logger): + parser = ElasticParserInterface() + + if self.parser is not None: + parser.reuse_parser(self.parser) + else: + self.parser = parser + + parser.parse(filepath, archive, logger) diff --git a/elasticparser/elastic_properties.py b/elasticparser/elastic_properties.py deleted file mode 100644 index 1f1a331..0000000 --- a/elasticparser/elastic_properties.py +++ /dev/null @@ -1,466 +0,0 @@ -import os -import pint -import numpy as np -from ase import Atoms - -from nomad.parsing.file_parser import Quantity, TextParser - - -class InfoParser(TextParser): - def __init__(self): - super().__init__(None) - - def init_quantities(self): - self._quantities = [ - Quantity( - 'order', r'\s*Order of elastic constants\s*=\s*([0-9]+)', repeats=False, - dtype=int), - Quantity( - 'calculation_method', r'\s*Method of calculation\s*=\s*([-a-zA-Z]+)\s*', - repeats=False), - Quantity( - 'code_name', r'\s*DFT code name\s*=\s*([-a-zA-Z]+)', repeats=False), - Quantity( - 'space_group_number', r'\s*Space-group number\s*=\s*([0-9]+)', repeats=False), - Quantity( - 'equilibrium_volume', r'\s*Volume of equilibrium unit cell\s*=\s*([0-9.]+)\s*', - unit='angstrom ** 3'), - Quantity( - 'max_strain', r'\s*Maximum Lagrangian strain\s*=\s*([0-9.]+)', repeats=False), - Quantity( - 'n_strains', r'\s*Number of distorted structures\s*=\s*([0-9]+)', repeats=False)] - - -class StructureParser(TextParser): - def __init__(self): - super().__init__(None) - - def init_quantities(self): - def get_sym_pos(val): - val = val.strip().replace('\n', '').split() - sym = [] - pos = [] - for i in range(0, len(val), 4): - sym.append(val[i + 3].strip()) - pos.append([float(val[j]) for j in range(i, i + 3)]) - sym_pos = dict(symbols=sym, positions=pos) - return sym_pos - - self._quantities = [ - Quantity( - 'cellpar', r'a\s*b\s*c\n([\d\.\s]+)\n\s*alpha\s*beta\s*gamma\n([\d\.\s]+)\n+', - repeats=False), - Quantity( - 'sym_pos', r'Atom positions:\n\n([\s\d\.A-Za-z]+)\n\n', - str_operation=get_sym_pos, repeats=False, convert=False)] - - -class DistortedParametersParser(TextParser): - def __init__(self): - super().__init__(None) - - def init_quantities(self): - self._quantities = [Quantity( - 'deformation', r'Lagrangian strain\s*=\s*\(([eta\s\d\.,]+)\)', - str_operation=lambda x: x.replace(',', '').split(), repeats=True, dtype=str)] - - -class FitParser(TextParser): - def __init__(self): - super().__init__(None) - - def init_quantities(self): - - def split_eta_val(val): - order, val = val.strip().split(' order fit.') - val = [float(v) for v in val.strip().split()] - return order, val[0::2], val[1::2] - - self._quantities = [Quantity( - 'fit', r'(\w+ order fit\.\n[\d.\s\neE\-\+]+)\n', repeats=True, convert=False, - str_operation=split_eta_val)] - - -class ElasticConstant2Parser(TextParser): - def __init__(self): - super().__init__(None) - - def init_quantities(self): - self._quantities = [ - Quantity( - 'voigt', r'Symmetry[\s\S]+\n\s*\n([C\d\s\n\(\)\-\+\/\*]+)\n', - shape=(6, 6), dtype=str, repeats=False), - Quantity( - 'elastic_constant', r'Elastic constant[\s\S]+in GPa\s*:\s*\n\n([\-\d\.\s\n]+)\n', - shape=(6, 6), dtype=float, unit='GPa', repeats=False), - Quantity( - 'compliance', r'Elastic compliance[\s\S]+in 1/GPa\s*:\s*\n\n([\-\d\.\s\n]+)\n', - shape=(6, 6), dtype=float, unit='1/GPa', repeats=False)] - - def str_to_modulus(val_in): - val_in = val_in.strip().split() - key = val_in[0] - unit = val_in[-1] if len(val_in) == 3 else None - val = float(val_in[1]) - val = pint.Quantity(val, unit) if unit is not None else val - return key, val - - self._quantities.append(Quantity( - 'modulus', r',\s*(\w+)\s*=\s*([\-\+\w\. ]+?)\n', str_operation=str_to_modulus, - repeats=True)) - - self._quantities.append(Quantity( - 'eigenvalues', - r'Eigenvalues of elastic constant \(stiffness\) matrix:\s*\n+([\-\d\.\n\s]+)\n', - unit='GPa', repeats=False)) - - -class ElasticConstant3Parser(TextParser): - def __init__(self): - super().__init__(None) - - def init_quantities(self): - def arrange_matrix(val): - val = val.strip().split('\n') - matrix = [v.strip().split() for v in val if v.strip()] - matrix = np.array(matrix).reshape((12, 18)) - arranged = [] - for i in range(2): - for j in range(3): - arranged.append( - matrix[i * 6: (i + 1) * 6, j * 6: (j + 1) * 6].tolist()) - return arranged - - self._quantities = [ - Quantity( - 'elastic_constant', r'\%\s*\n([\s0-6A-L]*)[\n\s\%1-6\-ij]*([\s0-6A-L]*)\n', - str_operation=arrange_matrix, dtype=str, repeats=False, convert=False), - Quantity( - 'cijk', r'(C\d\d\d)\s*=\s*([\-\d\.]+)\s*GPa', repeats=True, convert=False)] - - -class ElasticProperties: - def __init__(self): - self._mainfile = None - self.logger = None - self._deform_dirs = None - self._deform_dir_prefix = 'Dst' - self.info = InfoParser() - self.structure = StructureParser() - self.distorted_parameters = DistortedParametersParser() - self.fit = FitParser() - self.elastic_constant_2 = ElasticConstant2Parser() - self.elastic_constant_3 = ElasticConstant3Parser() - - @property - def mainfile(self): - return self._mainfile - - @mainfile.setter - def mainfile(self, val): - self._deform_dirs = None - self._mainfile = val - self.maindir = os.path.dirname(val) if val is not None else None - self.info.mainfile = val - - @property - def deformation_dirs(self): - if self._deform_dirs is None: - dirs = os.listdir(self.maindir) - self._deform_dirs = [ - os.path.join(self.maindir, d) for d in dirs if d.startswith(self._deform_dir_prefix)] - - return self._deform_dirs - - def get_references_to_calculations(self): - def output_file(dirname): - code = self.info.get('code_name', '').lower() - if code == 'exciting': - return os.path.join(dirname, 'INFO.OUT') - elif code == 'wien': - return os.path.join(dirname, '%s_Converged.scf' % os.path.basename(dirname)) - elif code == 'quantum': - return os.path.join(dirname, '%s.out' % os.path.basename(dirname)) - else: - return None - - references = [] - for deform_dir in self.deformation_dirs: - sub_dirs = os.listdir(deform_dir) - for sub_dir in sub_dirs: - calc_dir = os.path.join(deform_dir, sub_dir) - out_file = output_file(calc_dir) - if out_file is not None and os.path.isfile(out_file): - references.append(out_file) - - return references - - def get_structure_info(self): - path = os.path.join(self.maindir, 'sgroup.out') - if not os.path.isfile(path): - return - - self.structure.mainfile = path - - cellpar = self.structure.get('cellpar', None) - sym_pos = self.structure.get('sym_pos', {}) - - sym = sym_pos.get('symbols', None) - pos = sym_pos.get('positions', None) - - if cellpar is None or sym is None or pos is None: - return - - structure = Atoms(cell=cellpar, scaled_positions=pos, symbols=sym, pbc=True) - - positions = structure.get_positions() - positions = pint.Quantity(positions, 'angstrom') - cell = structure.get_cell() - cell = pint.Quantity(cell, 'angstrom') - - return sym, positions, cell - - def get_strain_energy(self): - strains, energies = [], [] - - for deform_dir in self.deformation_dirs: - filenames = [d for d in os.listdir(deform_dir) if d.endswith('Energy.dat')] - if not filenames: - continue - - path = os.path.join(deform_dir, filenames[-1]) - data = np.loadtxt(path).T - strains.append(list(data[0])) - # the peculiarity of the x_elastic_strain_diagram_values metainfo that it does - # not have the energy unit - energies.append(list(pint.Quantity(data[1], 'hartree').to('J').magnitude)) - - return strains, energies - - def get_strain_stress(self): - strains = {'Lagrangian-stress': [], 'Physical-stress': []} - stresses = {'Lagrangian-stress': [], 'Physical-stress': []} - - for deform_dir in self.deformation_dirs: - filenames = [d for d in os.listdir(deform_dir) if d.endswith('stress.dat')] - - for filename in filenames: - path = os.path.join(deform_dir, filename) - if not os.path.isfile(path): - continue - - with open(path) as f: - lines = f.readlines() - - strain, stress = [], [] - for line in lines: - val = line.strip().split() - if not val[0].strip().replace('.', '').isdecimal(): - continue - - strain.append(float(val[0])) - stress.append([float(v) for v in val[1:7]]) - - stype = filename.rstrip('.dat').split('_')[-1] - strains[stype].append(strain) - stresses[stype].append(stress) - - return strains, stresses - - def get_deformation_types(self): - path = os.path.join(self.maindir, 'Distorted_Parameters') - self.distorted_parameters.mainfile = path - return self.distorted_parameters.get('deformation') - - def _get_fit(self, path_dir, file_ext): - path_dir = os.path.join(self.maindir, path_dir) - - if not os.path.isdir(path_dir): - return - - paths = [p for p in os.listdir(path_dir) if p.endswith(file_ext)] - paths.sort() - - if not paths: - return - - eta, val = {}, {} - for path in paths: - self.fit.mainfile = os.path.join(path_dir, path) - fit_results = self.fit.get('fit', []) - for result in fit_results: - eta.setdefault(result[0], []) - val.setdefault(result[0], []) - eta[result[0]].append(result[1]) - val[result[0]].append(result[2]) - - return eta, val - - def get_energy_fit(self): - energy_fit = dict() - - for file_ext in ['d2E.dat', 'd3E.dat', 'ddE.dat']: - result = self._get_fit('Energy-vs-Strain', file_ext) - if result is None: - continue - - result = list(result) - result[1] = { - key: pint.Quantity( - val, 'GPa').to('Pa').magnitude for key, val in result[1].items()} - energy_fit['d2e'] = result - - result = self._get_fit('Energy-vs-Strain', 'CVe.dat') - if result is not None: - result = list(result) - result[1] = { - key: pint.Quantity( - val, 'hartree').to('J').magnitude for key, val in result[1].items()} - energy_fit['cross-validation'] = result - - return energy_fit - - def get_stress_fit(self): - stress_fit = dict() - stress_fit['dtn'] = [[]] * 6 - stress_fit['cross-validation'] = [[]] * 6 - - for strain_index in range(1, 7): - result = self._get_fit('Stress-vs-Strain', '%d_dS.dat' % strain_index) - if result is not None: - result[1] = {key: pint.Quantity(val, 'GPa') for key, val in result[1].items()} - stress_fit['dtn'][strain_index - 1] = result - - result = self._get_fit('Stress-vs-Strain', '%d_CVe.dat' % strain_index) - if result is not None: - result[1] = {key: pint.Quantity(val, 'hartree') for key, val in result[1].items()} - stress_fit['cross-validation'][strain_index - 1] = result - - return stress_fit - - def get_input(self): - paths = os.listdir(self.maindir) - path = None - order = self.info.get('order', 2) - for p in paths: - if 'ElaStic_' in p and p.endswith('.in') and str(order) in p: - path = p - break - - if path is None: - return - - calc_method = self.info.get('calculation_method') - - eta_ec = [] - fit_ec = [] - - def _is_number(var): - try: - float(var) - return True - except Exception: - return False - - path = os.path.join(self.maindir, path) - - with open(path) as f: - while True: - line = f.readline() - if not line: - break - - if calc_method.lower() == 'energy': - _, eta, fit = line.strip().split() - eta_ec.append(float(eta)) - fit_ec.append(int(fit)) - - elif calc_method.lower() == 'stress': - val = line.strip().split() - if not _is_number(val[0]): - eta_ec.append([float(val[i + 1]) for i in range(6)]) - else: - fit_ec.append([int(val[i]) for i in range(6)]) - - else: - pass - - return eta_ec, fit_ec - - def get_elastic_constants_order2(self): - path = os.path.join(self.maindir, 'ElaStic_2nd.out') - - self.elastic_constant_2.mainfile = path - - matrices = dict() - for key in ['voigt', 'elastic_constant', 'compliance']: - val = self.elastic_constant_2.get(key, None) - if val is not None: - matrices[key] = val - - moduli = dict() - for modulus in self.elastic_constant_2.get('modulus', []): - moduli[modulus[0]] = modulus[1] - - eigenvalues = self.elastic_constant_2.get('eigenvalues') - - return matrices, moduli, eigenvalues - - def get_elastic_constants_order3(self): - path = os.path.join(self.maindir, 'ElaStic_3rd.out') - self.elastic_constant_3.mainfile = path - - elastic_constant_str = self.elastic_constant_3.get('elastic_constant') - if elastic_constant_str is None: - return - - cijk = dict() - for element in self.elastic_constant_3.get('cijk', []): - cijk[element[0]] = float(element[1]) - - # formulas for the coefficients - coeff_A = cijk.get('C111', 0) + cijk.get('C112', 0) - cijk.get('C222', 0) - coeff_B = -(cijk.get('C115', 0) + 3 * cijk.get('C125', 0)) / 2 - coeff_C = (cijk.get('C114', 0) + 3 * cijk.get('C124', 0)) / 2 - coeff_D = -(2 * cijk.get('C111', 0) + cijk.get('C112', 0) - 3 * cijk.get('C222', 0)) / 4 - coeff_E = -cijk.get('C114', 0) - 2 * cijk.get('C124', 0) - coeff_F = -cijk.get('C115', 0) - 2 * cijk.get('C125', 0) - coeff_G = -(cijk.get('C115', 0) - cijk.get('C125', 0)) / 2 - coeff_H = (cijk.get('C114', 0) - cijk.get('C124', 0)) / 2 - coeff_I = (2 * cijk.get('C111', 0) - cijk.get('C112', 0) - cijk.get('C222', 0)) / 4 - coeff_J = (cijk.get('C113', 0) - cijk.get('C123', 0)) / 2 - coeff_K = -(cijk.get('C144', 0) - cijk.get('C155', 0)) / 2 - - space_group_number = self.info.get('space_group_number') - - if space_group_number <= 148: # rhombohedral II - coefficients = dict( - A=coeff_A, B=coeff_B, C=coeff_C, D=coeff_D, E=coeff_E, F=coeff_F, G=coeff_G, - H=coeff_H, I=coeff_I, J=coeff_J, K=coeff_K) - elif space_group_number <= 167: # rhombohedral I - coefficients = dict( - A=coeff_A, B=coeff_C, C=coeff_D, D=coeff_E, E=coeff_H, F=coeff_I, G=coeff_J, - H=coeff_K) - elif space_group_number <= 176: # hexagonal II - coefficients = dict( - A=coeff_A, B=coeff_D, C=coeff_I, D=coeff_J, E=coeff_K) - elif space_group_number <= 194: # hexagonal I - coefficients = dict( - A=coeff_A, B=coeff_D, C=coeff_I, D=coeff_J, E=coeff_K) - else: - coefficients = dict() - - # assign values to the elastic constant matrix from the independent ec and coefficients - elastic_constant = np.zeros((6, 6, 6)) - for i in range(6): - for j in range(6): - for k in range(6): - val = elastic_constant_str[i][j][k] - if val == '0': - elastic_constant[i][j][k] = 0 - elif val.isdigit(): - elastic_constant[i][j][k] = cijk['C%s' % val] - else: - elastic_constant[i][j][k] = coefficients.get(val, 0) - - return elastic_constant -- GitLab