From c0a152ed56d8f98687e0188382153cf23ba49a97 Mon Sep 17 00:00:00 2001 From: Alvin Noe Ladines <ladinesalvinnoe@gmail.com> Date: Thu, 3 Sep 2020 20:36:50 +0200 Subject: [PATCH] Switched to use nomad metainfo --- elasticparser/__init__.py | 19 +- elasticparser/__main__.py | 21 + elasticparser/elastic_parser.py | 177 ++ .../elastic_parser_input_exciting.py | 80 - elasticparser/elastic_parser_input_wien2k.py | 138 -- elasticparser/elastic_properties.py | 417 +++++ elasticparser/parser_elastic.py | 1453 ----------------- 7 files changed, 633 insertions(+), 1672 deletions(-) create mode 100644 elasticparser/__main__.py create mode 100644 elasticparser/elastic_parser.py delete mode 100644 elasticparser/elastic_parser_input_exciting.py delete mode 100644 elasticparser/elastic_parser_input_wien2k.py create mode 100644 elasticparser/elastic_properties.py delete mode 100644 elasticparser/parser_elastic.py diff --git a/elasticparser/__init__.py b/elasticparser/__init__.py index 9b38a23..018dcdb 100644 --- a/elasticparser/__init__.py +++ b/elasticparser/__init__.py @@ -12,4 +12,21 @@ # See the License for the specific language governing permissions and # limitations under the License. -from elasticparser.parser_elastic import ElasticParser +from .metainfo import m_env +from nomad.parsing.parser import FairdiParser +from elasticparser.elastic_parser import ElasticParserInterface + + +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')) + + def parse(self, filepath, archive, logger=None): + self._metainfo_env = m_env + + parser = ElasticParserInterface(filepath, archive, logger) + + parser.parse() diff --git a/elasticparser/__main__.py b/elasticparser/__main__.py new file mode 100644 index 0000000..47bbf8d --- /dev/null +++ b/elasticparser/__main__.py @@ -0,0 +1,21 @@ +# Copyright 2016-2018 Markus Scheidgen +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import sys +import json + +from elasticparser import ElasticParser + +if __name__ == "__main__": + archive = ElasticParser.main(sys.argv[1]) + json.dump(archive.m_to_dict(), sys.stdout, indent=2) \ No newline at end of file diff --git a/elasticparser/elastic_parser.py b/elasticparser/elastic_parser.py new file mode 100644 index 0000000..81a71a6 --- /dev/null +++ b/elasticparser/elastic_parser.py @@ -0,0 +1,177 @@ +import os +import numpy as np +import logging + +from nomadcore.unit_conversion import unit_conversion +from nomad.datamodel.metainfo.public import section_run, section_system,\ + section_single_configuration_calculation, section_method, section_calculation_to_calculation_refs,\ + Workflow, Elastic + +from elasticparser.metainfo.elastic import x_elastic_section_strain_diagrams,\ + x_elastic_section_fitting_parameters + +from elasticparser.elastic_properties import ElasticProperties + +a_to_m = unit_conversion.convert_unit_function('angstrom', 'meter') +ha_to_J = unit_conversion.convert_unit(1, 'hartree', 'J') +giga = 10 ** 9 + + +class ElasticParserInterface: + def __init__(self, filepath, archive, logger=None): + self.filepath = os.path.abspath(filepath) + self.archive = archive + self.logger = logger if logger is not None else logging + self.properties = ElasticProperties(self.filepath) + + def parse_strain(self): + sec_scc = self.archive.section_run[-1].section_single_configuration_calculation[-1] + method = self.properties.info['calculation_method'].lower() + + if method == 'energy': + strain, energy = self.properties.get_strain_energy() + n_strains = self.properties.info['n_strains'] + + sec_strain_diagram = sec_scc.m_create(x_elastic_section_strain_diagrams) + sec_strain_diagram.x_elastic_strain_diagram_type = 'energy' + sec_strain_diagram.x_elastic_strain_diagram_number_of_eta = len(strain[0]) + sec_strain_diagram.x_elastic_strain_diagram_eta_values = strain + sec_strain_diagram.x_elastic_strain_diagram_values = np.array(energy) * ha_to_J + + poly_fit_2 = int((n_strains - 1) / 2) + + poly_fit = { + '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() + for diagram_type in ['cross-validation', 'd2e']: + for fit_order in energy_fit[diagram_type][0].keys(): + sec_strain_diagram = sec_scc.m_create(x_elastic_section_strain_diagrams) + sec_strain_diagram.x_elastic_strain_diagram_type = diagram_type + sec_strain_diagram.x_elastic_strain_diagram_polynomial_fit_order = int(fit_order[:-2]) + sec_strain_diagram.x_elastic_strain_diagram_number_of_eta = poly_fit.get(fit_order, None) + sec_strain_diagram.x_elastic_strain_diagram_eta_values = energy_fit[diagram_type][0][fit_order] + convert = giga if diagram_type == 'd2e' else ha_to_J + sec_strain_diagram.x_elastic_strain_diagram_values = np.array(energy_fit[diagram_type][1][fit_order]) * convert + + elif method == 'stress': + strain, stress = self.properties.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)) + if not strain_i: + continue + + for si in range(6): + sec_strain_diagram = sec_scc.m_create(x_elastic_section_strain_diagrams) + sec_strain_diagram.x_elastic_strain_diagram_type = diagram_type + sec_strain_diagram.x_elastic_strain_diagram_stress_Voigt_component = si + 1 + sec_strain_diagram.x_elastic_strain_diagram_number_of_eta = len(strain_i[0]) + 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() + for diagram_type in ['cross-validation', 'dtn']: + for si in range(6): + for fit_order in stress_fit[diagram_type][si][0].keys(): + sec_strain_diagram = sec_scc.m_create(x_elastic_section_strain_diagrams) + sec_strain_diagram.x_elastic_strain_diagram_type = diagram_type + sec_strain_diagram.x_elastic_strain_diagram_stress_Voigt_component = si + 1 + sec_strain_diagram.x_elastic_strain_diagram_polynomial_fit_order = int(fit_order[:-2]) + sec_strain_diagram.x_elastic_strain_diagram_number_of_eta = poly_fit.get(fit_order, None) + sec_strain_diagram.x_elastic_strain_diagram_eta_values = stress_fit[diagram_type][si][0][fit_order] + convert = giga if diagram_type == 'dtn' else ha_to_J + sec_strain_diagram.x_elastic_strain_diagram_values = np.array(stress_fit[diagram_type][si][1][fit_order]) * convert + + def parse_elastic_constant(self): + sec_scc = self.archive.section_run[-1].section_single_configuration_calculation[-1] + + order = self.properties.info['order'] + + if order == 2: + matrices, moduli, eigenvalues = self.properties.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'] * giga + sec_scc.x_elastic_2nd_order_constants_compliance_matrix = matrices['compliance'] / giga + + sec_scc.x_elastic_Voigt_bulk_modulus = moduli.get('B_V', moduli.get('K_V')) * giga + sec_scc.x_elastic_Voigt_shear_modulus = moduli['G_V'] * giga + + sec_scc.x_elastic_Reuss_bulk_modulus = moduli.get('B_R', moduli.get('K_R')) * giga + sec_scc.x_elastic_Reuss_shear_modulus = moduli['G_R'] * giga + + sec_scc.x_elastic_Hill_bulk_modulus = moduli.get('B_H', moduli.get('K_H')) * giga + sec_scc.x_elastic_Hill_shear_modulus = moduli['G_H'] * giga + + sec_scc.x_elastic_Voigt_Young_modulus = moduli['E_V'] * giga + sec_scc.x_elastic_Voigt_Poisson_ratio = moduli['nu_V'] + sec_scc.x_elastic_Reuss_Young_modulus = moduli['E_R'] * giga + sec_scc.x_elastic_Reuss_Poisson_ratio = moduli['nu_R'] + sec_scc.x_elastic_Hill_Young_modulus = moduli['E_H'] * giga + sec_scc.x_elastic_Hill_Poisson_ratio = moduli['nu_H'] + + sec_scc.x_elastic_eigenvalues = eigenvalues * giga + + elif order == 3: + elastic_constant = self.properties.get_elastic_constants_order3() + + sec_scc.x_elastic_3rd_order_constants_matrix = elastic_constant + + def parse(self): + + sec_run = self.archive.m_create(section_run) + + sec_run.program_name = 'elastic' + sec_run.program_version = '1.0' + + sec_system = sec_run.m_create(section_system) + + symbols, positions, cell = self.properties.get_structure_info() + volume = self.properties.info['equilibrium_volume'] + volume = a_to_m(a_to_m(a_to_m(volume))) + + sec_system.atom_labels = symbols + sec_system.atom_positions = a_to_m(positions) + sec_system.simulation_cell = a_to_m(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_unit_cell_volume = volume + + sec_method = sec_run.m_create(section_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'] + + deformation_types = self.properties.get_deformation_types() + sec_method.x_elastic_deformation_types = deformation_types + sec_method.x_elastic_number_of_deformations = len(self.properties.deformation_dirs) + + references = self.properties.get_references_to_calculations() + sec_scc = sec_run.m_create(section_single_configuration_calculation) + for reference in references: + sec_calc_ref = sec_scc.m_create(section_calculation_to_calculation_refs) + sec_calc_ref.calculation_to_calculation_external_url = reference + sec_calc_ref.calculation_to_calculation_kind = 'source_calculation' + + fit_input = self.properties.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] + + self.parse_strain() + + self.parse_elastic_constant() + + sec_scc.single_configuration_to_calculation_method_ref = sec_method + sec_scc.single_configuration_calculation_to_system_ref = sec_system + + sec_workflow = self.archive.m_create(Workflow) + sec_workflow.workflow_type = 'elastic' + sec_elastic = sec_workflow.m_create(Elastic) + sec_elastic.calculation_method = self.properties.info['calculation_method'] + sec_elastic.elastic_constants_order = self.properties.info['order'] + sec_elastic.strain_maximum = self.properties.info['max_strain'] diff --git a/elasticparser/elastic_parser_input_exciting.py b/elasticparser/elastic_parser_input_exciting.py deleted file mode 100644 index 6ba57f6..0000000 --- a/elasticparser/elastic_parser_input_exciting.py +++ /dev/null @@ -1,80 +0,0 @@ -# Copyright 2017-2018 Lorenzo Pardini -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -import xml.sax -import logging -import numpy as np -from nomadcore.unit_conversion.unit_conversion import convert_unit_function -from nomadcore.unit_conversion.unit_conversion import convert_unit -from nomadcore.unit_conversion import unit_conversion - -class InputHandler(xml.sax.handler.ContentHandler): - - def __init__(self, backend): - self.backend = backend - self.inputSectionGIndex = -1 - self.basevect = [] - self.latticeDummy = '' - self.CurrentData = '' - self.atomCoor = [] - self.atomCoorDummy = [] - self.speciesfileDummy = '' - self.speciesfile = [] - self.scale = 1.0 - self.cell = [] - self.cellDummy = [] - - def endDocument(self): - bohr_to_m = convert_unit(1, "bohr", "m") - for i in range(0,len(self.cellDummy)): - for j in range(0,3): - self.cell[i].append(float(self.cellDummy[i][j])*self.scale*bohr_to_m) - self.backend.addValue("lattice_vectors", self.cell) - self.backend.addValue('atom_positions',self.atomCoor) - for i in range(0,len(self.atomCoor)): - self.speciesfile.append(self.speciesfileDummy) - self.backend.addValue("atom_labels", self.speciesfile) - def startElement(self, name, attrs): - self.CurrentData = name - if name == "crystal": - try: - self.scale = float(attrs.getValue('scale')) - except: - self.scale = 1.0 - elif name == 'species': - self.speciesfileDummy = attrs.getValue('speciesfile')[:-4] - elif name == 'atom': - self.atomCoorDummy = attrs.getValue('coord').split() - for j in range(0,3): - self.atomCoorDummy[j]=float(self.atomCoorDummy[j]) - self.atomCoor.append(self.atomCoorDummy) - else: - pass - - def endElement(self, name): - pass - - def characters(self, content): - if self.CurrentData == 'basevect': - self.latticeDummy = content - lattice = self.latticeDummy.split() - if lattice != []: - self.cellDummy.append(lattice) - self.cell.append([]) - else: - pass - -def parseInput(inF, backend): - handler = InputHandler(backend) - xml.sax.parse(inF, handler) diff --git a/elasticparser/elastic_parser_input_wien2k.py b/elasticparser/elastic_parser_input_wien2k.py deleted file mode 100644 index 107fa8f..0000000 --- a/elasticparser/elastic_parser_input_wien2k.py +++ /dev/null @@ -1,138 +0,0 @@ -from builtins import object -from nomadcore.simple_parser import mainFunction, CachingLevel -from nomadcore.simple_parser import SimpleMatcher as SM -from nomadcore.local_meta_info import loadJsonFile, InfoKindEl -import os, sys, json, logging -import numpy as np -import ase.geometry - - -################################################################ -# This is the subparser for the main WIEN2k input file (.struct) -################################################################ - -# Copyright 2016-2018 Daria M. Tomecka, Fawzi Mohamed -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -__author__ = "Daria M. Tomecka" -__maintainer__ = "Daria M. Tomecka" -__email__ = "tomeckadm@gmail.com;" -__date__ = "15/05/2017" - -########### Modified for elastic by Lorenzo Pardini ################## - -class Wien2kStructContext(object): - """context for wien2k struct parser""" - - def __init__(self): - self.parser = None - - def initialize_values(self): - """allows to reset values if the same superContext is used to parse different files""" - pass - - def startedParsing(self, path, parser): - """called when parsing starts""" - self.parser = parser - # allows to reset values if the same superContext is used to parse different files - self.initialize_values() - - def onClose_section_system(self, backend, gIndex, section): - # unit_cell - unit_cell_params = [] - for i in ['a', 'b', 'c']: - uci = section['x_elastic_wien2k_unit_cell_param_' + i] - #if uci is not None: - unit_cell_params.append(uci[0]) - for i in ['alfa', 'beta', 'gamma']: - uci = section['x_elastic_wien2k_angle_between_unit_axis_' + i] - # if uci is not None: - unit_cell_params.append(uci[0]) - - unit_cell = ase.geometry.cellpar_to_cell(unit_cell_params) - backend.addArrayValues('lattice_vectors', unit_cell) -# backend.addArrayValues("configuration_periodic_dimensions", np.ones(3, dtype=bool)) - - equiv_atoms = section["x_elastic_wien2k_section_equiv_atoms"] - #logging.error("section: %s", section) - labels = [] - pos = [] - - for eqAtoms in equiv_atoms: - label = eqAtoms["x_elastic_wien2k_atom_name"][0] - x = eqAtoms["x_elastic_wien2k_atom_pos_x"] - y = eqAtoms["x_elastic_wien2k_atom_pos_y"] - z = eqAtoms["x_elastic_wien2k_atom_pos_z"] - #logging.error("equiv_atoms: %s x %s y %s z %s",eqAtoms, x, y, z) - if len(x) != len(y) or len(x) != len(z): - raise Exception("incorrect parsing, different number of x,y,z components") - groupPos = [[x[i],y[i],z[i]] for i in range(len(x))] - nAt = len(groupPos) - labels += [label for i in range(nAt)] - pos += groupPos - backend.addValue("atom_labels", labels) - - backend.addArrayValues('atom_positions', np.asarray(pos)) - - - -# description of the input -def buildStructureMatchers(): - return SM( - name = 'root', - weak = True, - startReStr = "", - sections = ["section_run", "section_system"], - subMatchers = [ - SM(name = 'systemName', - startReStr = r"(?P<x_elastic_wien2k_system_nameIn>.*)"), - SM(r"\w+\s*LATTICE,NONEQUIV\.ATOMS.\s*(?P<x_elastic_wien2k_nonequiv_atoms>[0-9]+)"), - SM(r"(?P<x_elastic_wien2k_calc_mode>.*)"), - # SM(r"\s*(?P<x_wien2k_unit_cell_param_a>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*(?P<x_wien2k_unit_cell_param_b>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*(?P<x_wien2k_unit_cell_param_c>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*(?P<x_wien2k_angle_between_unit_axis_alfa>[-+0-9.eEdD]{9})\s*(?P<x_wien2k_angle_between_unit_axis_beta>[-+0-9.eEdD]{9})\s*(?P<x_wien2k_angle_between_unit_axis_gamma>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)"), - SM(r"\s*(?P<x_elastic_wien2k_unit_cell_param_a__angstrom>[-+0-9]*\.\d{0,6}){0,10}\s*(?P<x_elastic_wien2k_unit_cell_param_b__angstrom>[-+0-9]*\.\d{0,6}){0,10}\s*(?P<x_elastic_wien2k_unit_cell_param_c__angstrom>[-+0-9]*\.\d{0,6}){0,10}\s*(?P<x_elastic_wien2k_angle_between_unit_axis_alfa>[-+]?[0-9]*\.\d{0,6}){0,10}\s*(?P<x_elastic_wien2k_angle_between_unit_axis_beta>[-+]?[0-9]*\.\d{0,6}){0,10}\s*(?P<x_elastic_wien2k_angle_between_unit_axis_gamma>[-+]?[0-9]*\.\d*)"), - SM(r"\s*ATOM\s*[-0-9]+:\s*X=(?P<x_elastic_wien2k_atom_pos_x__angstrom>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*Y=(?P<x_elastic_wien2k_atom_pos_y__angstrom>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*Z=(?P<x_elastic_wien2k_atom_pos_z__angstrom>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)", - repeats=True, - sections=["x_elastic_wien2k_section_equiv_atoms"], - subMatchers=[ - SM(r"\s*[-0-9]+:\s*X=(?P<x_elastic_wien2k_atom_pos_x__angstrom>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*Y=(?P<x_elastic_wien2k_atom_pos_y__angstrom>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)\s*Z=(?P<x_elastic_wien2k_atom_pos_z__angstrom>[-+]?[0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?)", - repeats=True - ), - # SM(r"\s*(?P<atom>.{10})\s*NPT=\s*(?P<NPT>[0-9]+)\s*R0=(?P<r0>[0-9.]+)\s*RMT=\s*(?P<rmt>[0-9.]+)\s*Z:\s*(?P<z>[0-9.]+)",) - SM(r"\s*(?P<x_elastic_wien2k_atom_name>^.+)\s*NPT=\s*(?P<x_elastic_wien2k_NPT>[0-9]+)\s*R0=(?P<x_elastic_wien2k_R0>[0-9.]+)\s*RMT=\s*(?P<x_elastic_wien2k_RMT>[0-9.]+)\s*Z:\s*(?P<x_elastic_wien2k_atomic_number_Z>[0-9.]+)",) - ] - ) - ]) - -def get_cachingLevelForMetaName(metaInfoEnv, CachingLvl): - """Sets the caching level for the metadata. - - Args: - metaInfoEnv: metadata which is an object of the class InfoKindEnv in nomadcore.local_meta_info.py. - CachingLvl: Sets the CachingLevel for the sections k_band, run, and single_configuration_calculation. - This allows to run the parser without opening new sections. - - Returns: - Dictionary with metaname as key and caching level as value. - """ - # manually adjust caching of metadata - cachingLevelForMetaName = { - 'section_run': CachingLvl, - 'section_system': CachingLvl - } - cachingLevelForMetaName["x_elastic_wien2k_system_nameIn"] = CachingLevel.ForwardAndCache - cachingLevelForMetaName["x_elastic_wien2k_section_equiv_atoms"] = CachingLevel.ForwardAndCache - cachingLevelForMetaName["atom_labels"] = CachingLevel.ForwardAndCache - return cachingLevelForMetaName - -# loading metadata from nomad-meta-info/meta_info/nomad_meta_info/fhi_aims.nomadmetainfo.json diff --git a/elasticparser/elastic_properties.py b/elasticparser/elastic_properties.py new file mode 100644 index 0000000..214fb12 --- /dev/null +++ b/elasticparser/elastic_properties.py @@ -0,0 +1,417 @@ +import os +import numpy as np +from ase import Atoms + +from nomad.parsing.text_parser import Quantity, UnstructuredTextFileParser + + +class ElasticProperties: + def __init__(self, mainfile): + self._mainfile = mainfile + self.maindir = os.path.dirname(mainfile) + self._deform_dirs = None + self._deform_dir_prefix = 'Dst' + self._info = None + + @property + def mainfile(self): + return self._mainfile + + @mainfile.setter + def mainfile(self, val): + self._info = None + self._deform_dirs = None + self._mainfile = val + + @property + def info(self): + if self._info is None: + + quantities = [ + Quantity('order', r'\s*Order of elastic constants\s*=\s*([0-9]+)'), + Quantity('calculation_method', r'\s*Method of calculation\s*=\s*([-a-zA-Z]+)\s*'), + Quantity('code_name', r'\s*DFT code name\s*=\s*([-a-zA-Z]+)'), + Quantity('space_group_number', r'\s*Space-group number\s*=\s*([0-9]+)'), + Quantity('equilibrium_volume', r'\s*Volume of equilibrium unit cell\s*=\s*([-0-9.]+)\s*\[a.u\^3\]'), + Quantity('max_strain', r'\s*Maximum Lagrangian strain\s*=\s*([0-9.]+)'), + Quantity('n_strains', r'\s*Number of distorted structures\s*=\s*([0-9]+)') + ] + + parser = UnstructuredTextFileParser(self.mainfile, quantities) + + self._info = {} + for key, val in parser.items(): + self._info[key] = val[0] + + return self._info + + @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['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 + + 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 + + quantities = [ + Quantity('cellpar', r'a\s*b\s*c\n([\d\.\s]+)\n\s*alpha\s*beta\s*gamma\n([\d\.\s]+)\n+'), + Quantity('sym_pos', r'Atom positions:\n\n([\s\d\.A-Za-z]+)\n\n', str_operation=get_sym_pos) + ] + + parser = UnstructuredTextFileParser(path, quantities) + + cellpar = list(parser['cellpar'])[0] + sym_pos = list(parser['sym_pos'])[0] + + sym = list(sym_pos['symbols']) + pos = list(sym_pos['positions']) + + structure = Atoms(cell=cellpar, scaled_positions=pos, symbols=sym, pbc=True) + + return sym, structure.get_positions(), list(structure.get_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(data[0]) + energies.append(data[1]) + + 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') + if not os.path.isfile(path): + return + + quantities = [ + Quantity( + 'deformation', r'Lagrangian strain\s*=\s*\(([eta\s\d\.,]+)\)', + str_operation=lambda x: x.replace(',', '').split(), dtype=str)] + parser = UnstructuredTextFileParser(path, quantities) + + deformation_types = [list(d) for d in list(parser['deformation'])] + + return deformation_types + + 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 = {'2nd': [], '3rd': [], '4th': [], '5th': [], '6th': [], '7th': []} + val = {'2nd': [], '3rd': [], '4th': [], '5th': [], '6th': [], '7th': []} + + def split_eta_val(val): + val = val.strip().split() + return [val[0::2], val[1::2]] + + quantities = [ + Quantity( + order, r'%s order fit.\n([\d.\s\ne\-\+]+)\n' % order, str_operation=split_eta_val) + for order in eta.keys() + ] + + parser = UnstructuredTextFileParser(os.path.join(path_dir, paths[0]), quantities) + + for path in paths: + parser.mainfile = os.path.join(path_dir, path) + for key in eta.keys(): + if parser[key] is not None: + eta[key].append(list(parser[key][0][0])) + val[key].append(list(parser[key][0][1])) + + eta = {key: val for key, val in eta.items() if val} + val = {key: val for key, val in val.items() if val} + + 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 + + energy_fit['d2e'] = result + + result = self._get_fit('Energy-vs-Strain', 'CVe.dat') + if result is not None: + 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: + stress_fit['dtn'][strain_index - 1] = result + + result = self._get_fit('Stress-vs-Strain', '%d_CVe.dat' % strain_index) + if result is not None: + stress_fit['cross-validation'][strain_index - 1] = result + + return stress_fit + + def get_input(self): + paths = os.listdir(self.maindir) + path = None + for p in paths: + if 'ElaStic_' in p and p.endswith('.in'): + path = p + break + + if path is None: + return + + calc_method = self.info['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') + + def reshape(val): + return np.array(val.split()).reshape((6, 6)) + + quantities = [ + Quantity( + 'voigt', r'Symmetry of the second-order elastic constant matrix in Voigt notation.\s*[\s\S]+\n\n([C\d\s\n]+)\n', + str_operation=reshape, dtype=str), + Quantity( + 'elastic_constant', r'Elastic constant \(stiffness\) matrix in GPa:\s*\n\n([\-\d\.\s\n]+)\n', + str_operation=reshape, dtype=float), + Quantity( + 'compliance', r'Elastic compliance matrix in 1/GPa:\s*\n\n([\-\d\.\s\n]+)\n', + str_operation=reshape, dtype=float)] + + modulus_names = [ + 'B_V', 'K_V', 'G_V', 'B_R', 'K_R', 'G_R', 'B_H', 'K_H', 'G_H', 'E_V', + 'nu_V', 'E_R', 'nu_R', 'E_H', 'nu_H', 'AVR'] + + quantities += [ + Quantity(name, r'%s\s*=\s*([-+\d+\.]+)\s*' % name) for name in modulus_names] + + quantities += [ + Quantity( + 'eigenvalues', r'Eigenvalues of elastic constant \(stiffness\) matrix:\s*\n+([\-\d\.\n\s]+)\n')] + + parser = UnstructuredTextFileParser(path, quantities) + + matrices = dict( + voigt=list(parser['voigt'])[0], elastic_constant=list(parser['elastic_constant'])[0], + compliance=list(parser['compliance'])[0]) + + moduli = dict() + for name in modulus_names: + if parser[name] is not None: + moduli[name] = parser[name][0] + + eigenvalues = list(parser['eigenvalues'])[0] + + return matrices, moduli, eigenvalues + + def get_elastic_constants_order3(self): + path = os.path.join(self.maindir, 'ElaStic_3rd.out') + elastic_constant = [] + + 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 + + def get_cijk(val): + val = val.replace('=', '').replace('GPa', '') + cijk = dict() + for v in val.split('\n'): + v = v.split() + if len(v) == 2: + cijk[v[0].strip()] = float(v[1]) + return cijk + + space_group_number = self.info['space_group_number'] + + 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), + Quantity( + 'cijk', r'(C\d{3}\s*=\s*[C1-6\s=\n\-\.0-9GPa]*)\n\n', str_operation=get_cijk) + ] + + parsers = UnstructuredTextFileParser(path, quantities) + + elastic_constant_str = list(parsers['elastic_constant'])[0] + + cijk = list(parsers['cijk'])[0] + + # 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 + + 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 diff --git a/elasticparser/parser_elastic.py b/elasticparser/parser_elastic.py deleted file mode 100644 index 168f1ca..0000000 --- a/elasticparser/parser_elastic.py +++ /dev/null @@ -1,1453 +0,0 @@ -# Copyright 2017-2018 Lorenzo Pardini -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - -from builtins import object -import numpy as np -from nomadcore.unit_conversion.unit_conversion import convert_unit -from nomadcore.simple_parser import mainFunction, AncillaryParser, CachingLevel -from nomadcore.simple_parser import SimpleMatcher as SM -from nomadcore.local_meta_info import loadJsonFile, InfoKindEl -from nomadcore.unit_conversion import unit_conversion -import os, sys, json -import logging -import elasticparser.elastic_parser_input_exciting as elastic_parser_input_exciting -import elasticparser.elastic_parser_input_wien2k as elastic_parser_input_wien2k -from ase import Atoms -#from pathlib import Path - -def is_number(s): - try: - float(s) - return True - except ValueError: - return False - -class SampleContext(object): - - def __init__(self): -# self.mainFileUri = sys.argv[1] #exciting !!!!!!LOCAL HOME!!!!!!!! OKOKOKOK - self.mainFileUri = sys.argv[2] #exciting !!! FOR NOMAD URI nmd:// or sbt -> zip file!!!!!!!! OKOKKOOK - self.parser = None - self.mainFilePath = None - self.mainFile = None - self.etaEC = [] - self.fitEC = [] - self.SGN = 0 - self.secMethodIndex = None - self.secSystemIndex = None - self.sim_cell = [] - - def initialize_values(self): - """allows to reset values if the same superContext is used to parse different files""" - self.metaInfoEnv = self.parser.parserBuilder.metaInfoEnv - - def onOpen_section_system(self, backend, gIndex, section): - self.secSystemIndex = gIndex - - def startedParsing(self, path, parser): - """called when parsing starts""" - self.parser = parser - self.initialize_values() - - def onOpen_section_method(self, backend, gIndex, section): - if self.secMethodIndex is None: - self.secMethodIndex = gIndex - - def onClose_section_run(self, backend, gIndex, section): - backend.addValue('program_name', 'elastic') - backend.addValue('program_version', '1.0') - - def onClose_section_system(self, backend, gIndex, section): - backend.addArrayValues('configuration_periodic_dimensions', np.asarray([True, True, True])) - self.SGN = int(section["x_elastic_space_group_number"][0]) - mainFile = self.parser.fIn.fIn.name - dirPath = os.path.dirname(mainFile) #####exciting sbt -> zip file#### YES ?????????? sure???? check first OKOKOKKO - self.mainFile = self.parser.fIn.name - self.mainFilePath = self.mainFile[0:-12] -# dirPath = self.mainFileUri[0:-12] #####exciting LOCAL HOME or from NOMAD URI nmd:// ####### YES OKOKOKOK - for files in os.listdir(dirPath): - if files[-3:] == "xml": - inputFile = files - os.chdir(self.mainFilePath) - with open(inputFile) as f: - elastic_parser_input_exciting.parseInput(f, backend) - elif files[-6:] == "struct": - inputFile = files - os.chdir(self.mainFilePath) - structSuperContext = elastic_parser_input_wien2k.Wien2kStructContext() - structParser = AncillaryParser( - fileDescription = elastic_parser_input_wien2k.buildStructureMatchers(), - parser = self.parser, - cachingLevelForMetaName = elastic_parser_input_wien2k.get_cachingLevelForMetaName(self.metaInfoEnv, CachingLevel.PreOpenedIgnore), - superContext = structSuperContext) - - # ESPRESSO INPUT FILE TO BE ADDED - - with open(inputFile) as fIn: - structParser.parseFile(fIn) - - elif files[-3:] == ".in": ##### so far it works only for Rostam's calculations - if files != "ElaStic_2nd.in": - inputFile = files - atom_labels = [] - posX = [] - posY = [] - posZ = [] - coord = [] - lattice = [] - check = False - os.chdir(self.mainFilePath) - with open(inputFile) as g: - fromB = unit_conversion.convert_unit_function("bohr", "m") - while 1: - s = g.readline() - if not s: break - s = s.strip() - s = s.split() - if s[0] == "ibrav": ####### Rostam's: ibrav always 0 - ibrav = s[2] - elif s[0] == "celldm(1)": - alat = float(s[2]) - elif len(s) == 4: - atom_labels.append(s[0]) - posX.append(float(s[1])) - posY.append(float(s[2])) - posZ.append(float(s[3])) - elif len(s) == 3 and s[1] != "=": - if is_number(s[0]): - lattice.append([fromB(alat*float(s[0])),fromB(alat*float(s[1])),fromB(alat*float(s[2]))]) - else: - pass - - for i in range(len(atom_labels)): - coord.append([posX[i],posY[i],posZ[i]]) - cell = [[lattice[0][0],lattice[0][1],lattice[0][2]], - [lattice[1][0],lattice[1][1],lattice[1][2]], - [lattice[2][0],lattice[2][1],lattice[2][2]]] - self.sim_cell = cell - atoms = Atoms(atom_labels, coord, cell=[(1, 0, 0),(0, 1, 0),(0, 0, 1)]) - atoms.set_cell(self.sim_cell, scale_atoms=True) - coord = atoms.get_positions() - backend.addArrayValues('atom_positions', np.asarray(coord)) - backend.addArrayValues('atom_labels', np.asarray(atom_labels)) - backend.addValue("simulation_cell", cell) - - def onClose_section_method(self, backend, gIndex, section): - ha_per_joule = convert_unit(1, "hartree", "J") - giga = 10**9 - elCode = section['x_elastic_code'] - elasticGIndex = backend.openSection("section_single_configuration_calculation") - self.mainFilePath = self.mainFileUri[0:-12] - questa = os.getcwd() - mdr = float(section['x_elastic_max_lagrangian_strain'][0]) - ordr = int(section['x_elastic_elastic_constant_order'][0]) - nds = int(section['x_elastic_number_of_distorted_structures'][0]) - meth = section['x_elastic_calculation_method'][0] - polFit2 = (nds-1)/2 - polFit4 = polFit2 - 1 - polFit6 = polFit2 - 2 - polFit2Cross = polFit2 - 1 - polFit4Cross = polFit4 - 1 - polFit6Cross = polFit6 - 1 - polFit1 = (nds-1)/2 - polFit3 = polFit1 - 1 - polFit5 = polFit1 - 2 - polFit1Cross = polFit1 - 1 - polFit3Cross = polFit3 - 1 - polFit5Cross = polFit5 - 1 - ext_uri = [] - - i = 1 - while 1: - if (i<10): - Dstn = 'Dst0'+ str(i) - if (os.path.exists(Dstn) == True): - i += 1 - else: - break - else: - Dstn = 'Dst' + str(i) - if (os.path.exists(Dstn) == True): - i += 1 - else: - break - - defNum = i - 1 - ECs = defNum - - for j in range(1, ECs+1): - for i in range(1,nds+1): - if (j<10): - if (i<10): - if elCode[0] == 'exciting': - ext_uri.append(self.mainFilePath + 'Dst0' + str(j) + '/Dst0' + str(j) + '_0' + str(i) + '/INFO.OUT') - elif elCode[0] == 'WIEN': - ext_uri.append(self.mainFilePath + 'Dst0' + str(j) + '/Dst0' + str(j) + '_0' + str(i) + '/Dst0'+ str(j) + '_0' + str(i) + '_Converged.scf') - elif elCode[0] == 'QUANTUM': - ext_uri.append(self.mainFilePath + 'Dst0' + str(j) + '/Dst0' + str(j) + '_0' + str(i) + '/Dst0'+ str(j) + '_0' + str(i) + '.out') - else: - if elCode[0] == 'exciting': - ext_uri.append(self.mainFilePath + 'Dst0' + str(j) + '/Dst0' + str(j) + '_' + str(i) + '/INFO.OUT') - elif elCode[0] == 'WIEN': - ext_uri.append(self.mainFilePath + 'Dst0' + str(j) + '/Dst0' + str(j) + '_' + str(i) + '/Dst0' + str(j) + '_' + str(i) + '_Converged.scf') - elif elCode[0] == 'QUANTUM': - ext_uri.append(self.mainFilePath + 'Dst0' + str(j) + '/Dst0' + str(j) + '_' + str(i) + '/Dst0' + str(j) + '_' + str(i) + '.out') - else: - if (i<10): - if elCode[0] == 'exciting': - ext_uri.append(self.mainFilePath + 'Dst' + str(j) + '/Dst' + str(j) + '_0' + str(i) + '/INFO.OUT') - elif elCode[0] == 'WIEN': - ext_uri.append(self.mainFilePath + 'Dst' + str(j) + '/Dst' + str(j) + '_0' + str(i) + '/Dst' + str(j) + '_0' + str(i) + '_Converged.scf') - elif elCode[0] == 'QUANTUM': - ext_uri.append(self.mainFilePath + 'Dst' + str(j) + '/Dst' + str(j) + '_0' + str(i) + '/Dst' + str(j) + '_0' + str(i) + '.out') - else: - if elCode[0] == 'exciting': - ext_uri.append(self.mainFilePath + 'Dst' + str(j) + '/Dst' + str(j) + '_' + str(i) + '/INFO.OUT') - elif elCode[0] == 'WIEN': - ext_uri.append(self.mainFilePath + 'Dst' + str(j) + '/Dst' + str(j) + '_' + str(i) +'/Dst' + str(j) + '_' + str(i) + '_Converged.scf') - elif elCode[0] == 'QUANTUM': - ext_uri.append(self.mainFilePath + 'Dst' + str(j) + '/Dst' + str(j) + '_' + str(i) +'/Dst' + str(j) + '_' + str(i) + '.out') - for ref in ext_uri: - refGindex = backend.openSection("section_calculation_to_calculation_refs") - backend.addValue("calculation_to_calculation_external_url", ref) - backend.addValue("calculation_to_calculation_kind", "source_calculation") - backend.closeSection("section_calculation_to_calculation_refs", refGindex) - - energy = [] - eta = [] - LagrStress = [] - LagrStress_dummy = [] - physStress = [] - physStress_dummy = [] - for j in range(1, ECs+1): - if (j<10): - Dstn = 'Dst0'+ str(j) - eta.append([]) - energy.append([]) - LagrStress_dummy.append([]) - physStress_dummy.append([]) - else: - Dstn = 'Dst' + str(j) - eta.append([]) - energy.append([]) - LagrStress_dummy.append([]) - physStress_dummy.append([]) - - os.chdir(Dstn) - cur_dir = os.getcwd() - - if elCode[0] == 'exciting': - try: - f = open(Dstn+'-Energy.dat', 'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - dummy_eta, dummy_energy = s.split() - eta[-1].append(float(dummy_eta)) - energy[-1].append(float(dummy_energy)*ha_per_joule) - os.chdir('../') - except: - pass - try: - f = open(Dstn+'_Energy.dat', 'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - dummy_eta, dummy_energy = s.split() - eta[-1].append(float(dummy_eta)) - energy[-1].append(float(dummy_energy)*ha_per_joule) - os.chdir('../') - except: - pass - - elif elCode[0] == 'WIEN': - f = open(Dstn+'_Energy.dat', 'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - dummy_eta, dummy_energy = s.split() - eta[-1].append(float(dummy_eta)) - energy[-1].append(float(dummy_energy)*ha_per_joule) - os.chdir('../') - - elif (elCode[0] == 'QUANTUM' or elCode[0] == 'Quantum') and meth == 'Energy': - f = open(Dstn+'_Energy.dat', 'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - dummy_eta, dummy_energy = s.split() - eta[-1].append(float(dummy_eta)) - energy[-1].append(float(dummy_energy)*ha_per_joule) - os.chdir('../') - - elif elCode[0] == 'QUANTUM' and meth == 'Stress': - f = open(Dstn+'_Lagrangian-stress.dat', 'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - s = s.split() - if is_number(s[0]): - dummy_eta = s[0] - dummy_LS1 = s[1] - dummy_LS2 = s[2] - dummy_LS3 = s[3] - dummy_LS4 = s[4] - dummy_LS5 = s[5] - dummy_LS6 = s[6] - eta[-1].append(float(dummy_eta)) - LagrStress_dummy[-1].append([float(dummy_LS1),float(dummy_LS2),float(dummy_LS3),float(dummy_LS4),float(dummy_LS5),float(dummy_LS6)]) - - g = open(Dstn+'_Physical-stress.dat', 'r') - while 1: - s = g.readline() - if not s: break - s = s.strip() - s = s.split() - if is_number(s[0]): - dummy_PS1 = s[1] - dummy_PS2 = s[2] - dummy_PS3 = s[3] - dummy_PS4 = s[4] - dummy_PS5 = s[5] - dummy_PS6 = s[6] - physStress_dummy[-1].append([float(dummy_PS1),float(dummy_PS2),float(dummy_PS3),float(dummy_PS4),float(dummy_PS5),float(dummy_PS6)]) - os.chdir('../') - - else: - os.chdir('../') - - defTyp = [] - - f = open('Distorted_Parameters','r') - - while 1: - s = f.readline() - if not s: break - s = s.strip() - if 'Lagrangian' in s: - defTyp.append([]) - s = s.split("(") - s = s[-1].split(")") - s = s[0].split(",") - for i in range(0,6): - s[i] = s[i].strip() - defTyp[-1].append(s[i]) - - f.close() - prova = os.listdir('.') - if 'Energy-vs-Strain' in prova: - os.chdir('Energy-vs-Strain') - - d2e6_val = [] - d2e4_val = [] - d2e2_val = [] - d2e6_eta = [] - d2e4_eta = [] - d2e2_eta = [] - d2e_val_tot = [] - d2e_eta_tot = [] - - for i in range(1, ECs+1): - d2e_val_tot.append([]) - d2e_eta_tot.append([]) - if (i<10): - if ordr == 2: - Dstna = 'Dst0'+ str(i) + '_d2E.dat' - Dstnb = 'Dst0'+ str(i) + '_ddE.dat' - elif ordr == 3: - Dstna = 'Dst0'+ str(i) + '_d3E.dat' -# Dstnb = 'Dst0'+ str(i) + '-d3E.dat' - try: - f = open (Dstna,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - d2e_val_tot[-1].append([]) - d2e_eta_tot[-1].append([]) - elif len(s) >= 30: - d2e_eta, d2e_values = s.split() - d2e_val_tot[-1][-1].append(float(d2e_values)*giga) - d2e_eta_tot[-1][-1].append(float(d2e_eta)) - f.close() - except: - pass - try: - f = open (Dstnb,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - d2e_val_tot[-1].append([]) - d2e_eta_tot[-1].append([]) - elif len(s) >= 30: - d2e_eta, d2e_values = s.split() - d2e_val_tot[-1][-1].append(float(d2e_values)*giga) - d2e_eta_tot[-1][-1].append(float(d2e_eta)) - f.close() - except: - pass - else: - if ordr == 2: - Dstna = 'Dst' + str(i) + '_d2E.dat' - Dstnb = 'Dst' + str(i) + '_ddE.dat' - elif ordr == 3: - Dstna = 'Dst'+ str(i) + '_d3E.dat' - try: - f = open (Dstna,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - d2e_val_tot[-1].append([]) - d2e_eta_tot[-1].append([]) - elif len(s) >= 30: - d2e_eta, d2e_values = s.split() - d2e_val_tot[-1][-1].append(float(d2e_values)*giga) - d2e_eta_tot[-1][-1].append(float(d2e_eta)) - f.close() - except: - pass - try: - f = open (Dstnb,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - d2e_val_tot[-1].append([]) - d2e_eta_tot[-1].append([]) - elif len(s) >= 30: - d2e_eta, d2e_values = s.split() - d2e_val_tot[-1][-1].append(float(d2e_values)*giga) - d2e_eta_tot[-1][-1].append(float(d2e_eta)) - f.close() - except: - pass - - d2e6_val.append(d2e_val_tot[i-1][0]) - d2e4_val.append(d2e_val_tot[i-1][1]) - d2e2_val.append(d2e_val_tot[i-1][2]) - d2e6_eta.append(d2e_eta_tot[i-1][0]) - d2e4_eta.append(d2e_eta_tot[i-1][1]) - d2e2_eta.append(d2e_eta_tot[i-1][2]) - CrossVal6_val = [] - CrossVal4_val = [] - CrossVal2_val = [] - CrossVal_val_tot = [] - - CrossVal6_eta = [] - CrossVal4_eta = [] - CrossVal2_eta = [] - CrossVal_eta_tot = [] - - for i in range(1, ECs+1): - CrossVal_val_tot.append([]) - CrossVal_eta_tot.append([]) - if (i<10): - DstnCV = 'Dst0'+ str(i) + '_CVe.dat' - f = open (DstnCV,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - CrossVal_val_tot[-1].append([]) - CrossVal_eta_tot[-1].append([]) - elif len(s) >= 20 and s.split()[0] != '#': - CrossVal_eta, CrossVal_values = s.split() - CrossVal_val_tot[-1][-1].append(float(CrossVal_values)*ha_per_joule) - CrossVal_eta_tot[-1][-1].append(float(CrossVal_eta)) - f.close() - else: - DstnCV = 'Dst' + str(i) + '_CVe.dat' - f = open (Dstn,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - CrossVal_val_tot[-1].append([]) - CrossVal_eta_tot[-1].append([]) - elif len(s) >= 20 and s.split()[0] != '#': - CrossVal_eta, CrossVal_values = s.split() - CrossVal_val_tot[-1][-1].append(float(CrossVal_values)*ha_per_joule) - CrossVal_eta_tot[-1][-1].append(float(CrossVal_eta)) - f.close() - CrossVal6_val.append(CrossVal_val_tot[i-1][0]) - CrossVal4_val.append(CrossVal_val_tot[i-1][1]) - CrossVal2_val.append(CrossVal_val_tot[i-1][2]) - CrossVal6_eta.append(CrossVal_eta_tot[i-1][0]) - CrossVal4_eta.append(CrossVal_eta_tot[i-1][1]) - CrossVal2_eta.append(CrossVal_eta_tot[i-1][2]) - - os.chdir('../') - - elif 'Stress-vs-Strain' in prova: - - os.chdir('Stress-vs-Strain') - - dS5_val = [[],[],[],[],[],[]] - dS3_val = [[],[],[],[],[],[]] - dS1_val = [[],[],[],[],[],[]] - dS5_eta = [[],[],[],[],[],[]] - dS3_eta = [[],[],[],[],[],[]] - dS1_eta = [[],[],[],[],[],[]] - dS_val_tot = [[],[],[],[],[],[]] - dS_eta_tot = [[],[],[],[],[],[]] - string = [] - stringCV = [] - - for i in range(6): - string.append('Dstn'+str(i+1)) - stringCV.append('DstnCV'+str(i+1)) - - for Dstn in string: - j = string.index(Dstn) - for i in range(1, ECs+1): - dS_val_tot[string.index(Dstn)].append([]) - dS_eta_tot[string.index(Dstn)].append([]) - if (i<10): - Dstn = 'Dst0'+ str(i) + '_LS' + str(string.index(Dstn)+1) + '_dS.dat' - f = open (Dstn,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - dS_val_tot[j][-1].append([]) - dS_eta_tot[j][-1].append([]) - elif len(s) >= 30: - dS_eta, dS_values = s.split() - dS_val_tot[j][-1][-1].append(float(dS_values)*giga) - dS_eta_tot[j][-1][-1].append(float(dS_eta)) - f.close() - else: - Dstn = 'Dst'+ str(i) + '_LS' + str(string.index(Dstn)+1) + '_dS.dat' - f = open (Dstn,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - dS_val_tot[j][-1].append([]) - dS_eta_tot[j][-1].append([]) - elif len(s) >= 30: - dS_eta, dS_values = s.split() - dS_val_tot[j][-1][-1].append(float(dS_values)*giga) - dS_eta_tot[j][-1][-1].append(float(dS_eta)) - f.close() - - dS5_val[j].append(dS_val_tot[j][i-1][0]) - dS3_val[j].append(dS_val_tot[j][i-1][1]) - dS1_val[j].append(dS_val_tot[j][i-1][2]) - dS5_eta[j].append(dS_eta_tot[j][i-1][0]) - dS3_eta[j].append(dS_eta_tot[j][i-1][1]) - dS1_eta[j].append(dS_eta_tot[j][i-1][2]) - - CrossVal5_val = [[],[],[],[],[],[]] - CrossVal3_val = [[],[],[],[],[],[]] - CrossVal1_val = [[],[],[],[],[],[]] - CrossVal_val_tot = [[],[],[],[],[],[]] - - CrossVal5_eta = [[],[],[],[],[],[]] - CrossVal3_eta = [[],[],[],[],[],[]] - CrossVal1_eta = [[],[],[],[],[],[]] - CrossVal_eta_tot = [[],[],[],[],[],[]] - - for DstnCV in stringCV: - j = stringCV.index(DstnCV) - for i in range(1, ECs+1): - CrossVal_val_tot[stringCV.index(DstnCV)].append([]) - CrossVal_eta_tot[stringCV.index(DstnCV)].append([]) - if (i<10): - DstnCV = 'Dst0'+ str(i) + '_LS' + str(stringCV.index(DstnCV)+1) + '_CVe.dat' - f = open (DstnCV,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - CrossVal_val_tot[j][-1].append([]) - CrossVal_eta_tot[j][-1].append([]) - elif len(s) >= 20 and s.split()[0] != '#': - CrossVal_eta, CrossVal_values = s.split() - CrossVal_val_tot[j][-1][-1].append(float(CrossVal_values)*ha_per_joule) - CrossVal_eta_tot[j][-1][-1].append(float(CrossVal_eta)) - f.close() - else: - DstnCV = 'Dst'+ str(i) + '_LS' + str(stringCV.index(DstnCV)+1) + '_CVe.dat' - f = open (DstnCV,'r') - while 1: - s = f.readline() - if not s: break - s = s.strip() - if "order" in s.split(): - CrossVal_val_tot[j][-1].append([]) - CrossVal_eta_tot[j][-1].append([]) - elif len(s) >= 20 and s.split()[0] != '#': - CrossVal_eta, CrossVal_values = s.split() - CrossVal_val_tot[j][-1][-1].append(float(CrossVal_values)*ha_per_joule) - CrossVal_eta_tot[j][-1][-1].append(float(CrossVal_eta)) - f.close() - CrossVal5_val[j].append(CrossVal_val_tot[j][i-1][0]) - CrossVal3_val[j].append(CrossVal_val_tot[j][i-1][1]) - CrossVal1_val[j].append(CrossVal_val_tot[j][i-1][2]) - CrossVal5_eta[j].append(CrossVal_eta_tot[j][i-1][0]) - CrossVal3_eta[j].append(CrossVal_eta_tot[j][i-1][1]) - CrossVal1_eta[j].append(CrossVal_eta_tot[j][i-1][2]) - os.chdir('../') - - else: - pass - - if ordr == 2: - f = open ('ElaStic_'+str(ordr)+'nd.in','r') - elif ordr == 3: - f = open ('ElaStic_'+str(ordr)+'rd.in','r') - - EC_eigen = [] - - if meth == 'Energy': - for i in range(1, ECs+1): - s = f.readline() - s = s.strip() - dummy, etaEC_dummy, fitEC_dummy = s.split() - self.etaEC.append(float(etaEC_dummy)) - self.fitEC.append(int(fitEC_dummy)) - - elif meth == 'Stress': - while 1: - s = f.readline() - if not s: break - s = s.strip() - s = s.split() - if not is_number(s[0]): - self.etaEC.append([]) - for i in range(6): self.etaEC[-1].append(float(s[i+1])) - else: - self.fitEC.append([]) - for i in range(6): self.fitEC[-1].append(int(s[i])) - else: - pass - - f.close() - - if ordr == 2: - f = open ('ElaStic_'+str(ordr)+'nd.out','r') - - allMat = [[],[],[],[],[],[]] - voigtMat = [[],[],[],[],[],[]] - ECMat = [[],[],[],[],[],[]] - complMat = [[],[],[],[],[],[]] - - while 1: - s = f.readline() - if not s: break - s = s.strip() - s = s.split() - if len(s) == 1: - try: float(s[0]) - except ValueError: - continue - else: - EC_eigen.append(float(s[0])*giga) - elif "B_V" in s: - B_V = float(s[5])*giga - elif "K_V" in s: - B_V = float(s[5])*giga - elif "G_V" in s: - G_V = float(s[5])*giga - elif "B_R" in s: - B_R = float(s[5])*giga - elif "K_R" in s: - B_R = float(s[5])*giga - elif "G_R" in s: - G_R = float(s[5])*giga - elif "B_H" in s: - B_H = float(s[5])*giga - elif "K_H" in s: - B_H = float(s[5])*giga - elif "G_H" in s: - G_H = float(s[5])*giga - elif "E_V" in s: - E_V = float(s[5])*giga - elif "nu_V" in s: - nu_V = float(s[5]) - elif "E_R" in s: - E_R = float(s[5])*giga - elif "nu_R" in s: - nu_R = float(s[5]) - elif "E_H" in s: - E_H = float(s[5])*giga - elif "nu_H" in s: - nu_H = float(s[5]) - elif len(s) == 6 and s[0] != "Elastic" and s[0] != "Eigenvalues": - for i in range(0,6): - allMat[i].append(s[i]) - elif "AVR" in s: - AVR = float(s[6]) - - f.close() - - for i in range(0,6): - voigtMat[i] = allMat[i][0:6] - ECMat[i] = allMat[i][6:12] - complMat[i] = allMat[i][12:18] - - for i in range(0,6): - for j in range(0,6): - voigtMat[i][j] = voigtMat[j][i] - ECMat[i][j] = float(ECMat[j][i])*giga - complMat[i][j] = float(complMat[j][i])/giga - - if meth == 'Energy': - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "energy") - backend.addValue("x_elastic_strain_diagram_number_of_eta", len(eta[0])) - backend.addValue("x_elastic_strain_diagram_eta_values", eta) - backend.addValue("x_elastic_strain_diagram_values", energy) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "cross-validation") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 2) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit2Cross) - backend.addValue("x_elastic_strain_diagram_eta_values", CrossVal2_eta) - backend.addValue("x_elastic_strain_diagram_values", CrossVal2_val) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "cross-validation") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 4) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit4Cross) - backend.addValue("x_elastic_strain_diagram_eta_values", CrossVal4_eta) - backend.addValue("x_elastic_strain_diagram_values", CrossVal4_val) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "cross-validation") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 6) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit6Cross) - backend.addValue("x_elastic_strain_diagram_eta_values", CrossVal6_eta) - backend.addValue("x_elastic_strain_diagram_values", CrossVal6_val) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "d2e") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 2) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit2) - backend.addValue("x_elastic_strain_diagram_eta_values", d2e2_eta) - backend.addValue("x_elastic_strain_diagram_values", d2e2_val) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "d2e") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 4) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit4) - backend.addValue("x_elastic_strain_diagram_eta_values", d2e4_eta) - backend.addValue("x_elastic_strain_diagram_values", d2e4_val) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "d2e") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 6) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit6) - backend.addValue("x_elastic_strain_diagram_eta_values", d2e6_eta) - backend.addValue("x_elastic_strain_diagram_values", d2e6_val) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elif meth == 'Stress': - for k in range(6): - LagrStress.append([]) - physStress.append([]) - for j in range(len(LagrStress_dummy)): - LagrStress[-1].append([]) - physStress[-1].append([]) - for i in range(len(LagrStress_dummy[j])): - LagrStress[k][j].append(LagrStress_dummy[j][i][k]) - physStress[k][j].append(physStress_dummy[j][i][k]) - - for i in range(0,6): - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "Lagrangian-stress") - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_number_of_eta", len(eta[0])) - backend.addValue("x_elastic_strain_diagram_eta_values", eta) - backend.addValue("x_elastic_strain_diagram_values", LagrStress[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "Physical-stress") - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_number_of_eta", len(eta[0])) - backend.addValue("x_elastic_strain_diagram_eta_values", eta) - backend.addValue("x_elastic_strain_diagram_values", physStress[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "cross-validation") - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 1) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit1Cross) - backend.addValue("x_elastic_strain_diagram_eta_values", CrossVal1_eta[i]) - backend.addValue("x_elastic_strain_diagram_values", CrossVal1_val[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "cross-validation") - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 3) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit3Cross) - backend.addValue("x_elastic_strain_diagram_eta_values", CrossVal3_eta[i]) - backend.addValue("x_elastic_strain_diagram_values", CrossVal3_val[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "cross-validation") - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 5) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit5Cross) - backend.addValue("x_elastic_strain_diagram_eta_values", CrossVal5_eta[i]) - backend.addValue("x_elastic_strain_diagram_values", CrossVal5_val[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "dtn") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 1) - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit1) - backend.addValue("x_elastic_strain_diagram_eta_values", dS1_eta[i]) - backend.addValue("x_elastic_strain_diagram_values", dS1_val[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "dtn") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 3) - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit3) - backend.addValue("x_elastic_strain_diagram_eta_values", dS3_eta[i]) - backend.addValue("x_elastic_strain_diagram_values", dS3_val[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "dtn") - backend.addValue("x_elastic_strain_diagram_polynomial_fit_order", 5) - backend.addValue("x_elastic_strain_diagram_stress_Voigt_component", int(i+1)) - backend.addValue("x_elastic_strain_diagram_number_of_eta", polFit5) - backend.addValue("x_elastic_strain_diagram_eta_values", dS5_eta[i]) - backend.addValue("x_elastic_strain_diagram_values", dS5_val[i]) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - - else: - pass - - backend.addValue('x_elastic_2nd_order_constants_notation_matrix',voigtMat) - backend.addValue('x_elastic_2nd_order_constants_matrix',ECMat) - backend.addValue('x_elastic_2nd_order_constants_compliance_matrix',complMat) - backend.addValue('x_elastic_Voigt_bulk_modulus',B_V) - backend.addValue('x_elastic_Voigt_shear_modulus',G_V) - backend.addValue('x_elastic_Reuss_bulk_modulus',B_R) - backend.addValue('x_elastic_Reuss_shear_modulus',G_R) - backend.addValue('x_elastic_Hill_bulk_modulus',B_H) - backend.addValue('x_elastic_Hill_shear_modulus',G_H) - backend.addValue('x_elastic_Voigt_Young_modulus',E_V) - backend.addValue('x_elastic_Voigt_Poisson_ratio',nu_V) - backend.addValue('x_elastic_Reuss_Young_modulus',E_R) - backend.addValue('x_elastic_Reuss_Poisson_ratio',nu_R) - backend.addValue('x_elastic_Hill_Young_modulus',E_H) - backend.addValue('x_elastic_Hill_Poisson_ratio',nu_H) - backend.addValue('x_elastic_eigenvalues',EC_eigen) - backend.closeSection("section_single_configuration_calculation", elasticGIndex) - backend.addValue("x_elastic_deformation_types", defTyp) - backend.addValue("x_elastic_number_of_deformations", defNum) - if meth == "Energy": - elasticPIndex = backend.openSection("x_elastic_section_fitting_parameters") - backend.addValue("x_elastic_fitting_parameters_eta", self.etaEC) - backend.addValue("x_elastic_fitting_parameters_polynomial_order", self.fitEC) - backend.closeSection("x_elastic_section_fitting_parameters", elasticPIndex) - else: - pass - - elif ordr == 3: - f = open ('ElaStic_'+str(ordr)+'rd.out','r') - - ECmat = [] - for i in range(0,6): - ECmat.append([]) - for j in range(0,6): - ECmat[i].append([]) - for k in range(0,6): - ECmat[i][j].append([]) - ECmat[i][j][k] = int(0) - - while 1: - s = f.readline() - if not s: break - s = s.strip() - s = s.split() - if len(s) == 4: - if s[0] == 'C111': - C111 = float(s[2]) - elif s[0] == 'C112': - C112 = float(s[2]) - elif s[0] == 'C123': - C123 = float(s[2]) - elif s[0] == 'C144': - C144 = float(s[2]) - elif s[0] == 'C155': - C155 = float(s[2]) - elif s[0] == 'C456': - C456 = float(s[2]) - elif s[0] == 'C113': - C113 = float(s[2]) - elif s[0] == 'C166': - C166 = float(s[2]) - elif s[0] == 'C133': - C133 = float(s[2]) - elif s[0] == 'C222': - C222 = float(s[2]) - elif s[0] == 'C333': - C333 = float(s[2]) - elif s[0] == 'C344': - C344 = float(s[2]) - elif s[0] == 'C116': - C116 = float(s[2]) - elif s[0] == 'C145': - C145 = float(s[2]) - elif s[0] == 'C114': - C114 = float(s[2]) - elif s[0] == 'C124': - C124 = float(s[2]) - elif s[0] == 'C134': - C134 = float(s[2]) - elif s[0] == 'C444': - C444 = float(s[2]) - - - if(149 <= self.SGN and self.SGN <= 167): # Rhombohedral I - LC = 'RI' - ECs= 14 - - A = C111-C222+C112 - B = float(3/4)*C222-float(1/2)*C111-float(1/4)*C112 - C = float(1/2)*C111-float(1/4)*C222-float(1/4)*C112 - D = float(1/2)*(C113-C123) - E = float(1/2)*(C155-C144) - F = float(1/2)*(C114+float(3)*C124) - G =-C114-float(2)*C124 - H = float(1/2)*(C114-C124) - - ECmat[0][0][0] = C111 - ECmat[0][0][1] = C112 - ECmat[0][0][2] = C113 - ECmat[0][0][3] = C114 - ECmat[0][1][0] = C112 - ECmat[0][1][1] = A - ECmat[0][1][2] = C123 - ECmat[0][1][3] = C124 - ECmat[0][2][0] = C113 - ECmat[0][2][1] = C123 - ECmat[0][2][2] = C133 - ECmat[0][2][3] = C134 - ECmat[0][3][0] = C114 - ECmat[0][3][1] = C124 - ECmat[0][3][2] = C134 - ECmat[0][3][3] = C144 - ECmat[0][4][4] = C155 - ECmat[0][4][5] = F - ECmat[0][5][4] = F - ECmat[0][5][5] = B - - ECmat[1][0][0] = C112 - ECmat[1][0][1] = A - ECmat[1][0][2] = C123 - ECmat[1][0][3] = C124 - ECmat[1][1][0] = A - ECmat[1][1][1] = C222 - ECmat[1][1][2] = C113 - ECmat[1][1][3] = G - ECmat[1][2][0] = C123 - ECmat[1][2][1] = C113 - ECmat[1][2][2] = C133 - ECmat[1][2][3] = -C134 - ECmat[1][3][0] = C124 - ECmat[1][3][1] = G - ECmat[1][3][2] = -C134 - ECmat[1][3][3] = C155 - ECmat[1][4][4] = C144 - ECmat[1][4][5] = H - ECmat[1][5][4] = H - ECmat[1][5][5] = C - - ECmat[2][0][0] = C113 - ECmat[2][0][1] = C123 - ECmat[2][0][2] = C133 - ECmat[2][0][3] = C134 - ECmat[2][1][0] = C123 - ECmat[2][1][1] = C113 - ECmat[2][1][2] = C133 - ECmat[2][1][3] = -C134 - ECmat[2][2][0] = C133 - ECmat[2][2][1] = C133 - ECmat[2][2][2] = C133 - ECmat[2][3][0] = C134 - ECmat[2][3][1] = -C134 - ECmat[2][3][3] = C344 - ECmat[2][4][4] = C344 - ECmat[2][4][5] = C134 - ECmat[2][5][4] = C134 - ECmat[2][5][5] = D - - ECmat[3][0][0] = C114 - ECmat[3][0][1] = C124 - ECmat[3][0][2] = C134 - ECmat[3][0][3] = C144 - ECmat[3][1][0] = C12 - ECmat[3][1][1] = G - ECmat[3][1][2] = -C134 - ECmat[3][1][3] = C155 - ECmat[3][2][0] = C134 - ECmat[3][2][1] = -C134 - ECmat[3][2][3] = C344 - ECmat[3][3][0] = C144 - ECmat[3][3][1] = C155 - ECmat[3][3][2] = C344 - ECmat[3][3][3] = C444 - ECmat[3][4][4] = -C444 - ECmat[3][4][5] = E - ECmat[3][5][4] = E - ECmat[3][5][5] = C124 - - ECmat[4][0][4] = C155 - ECmat[4][0][5] = F - ECmat[4][1][4] = C144 - ECmat[4][1][5] = H - ECmat[4][2][4] = C344 - ECmat[4][2][5] = C134 - ECmat[4][3][5] = E - ECmat[4][4][0] = C155 - ECmat[4][4][1] = C144 - ECmat[4][4][2] = C344 - ECmat[4][5][0] = F - ECmat[4][5][1] = H - ECmat[4][5][2] = C134 - ECmat[4][5][3] = E - - ECmat[5][0][4] = F - ECmat[5][0][5] = B - ECmat[5][1][4] = H - ECmat[5][1][5] = C - ECmat[5][2][4] = C134 - ECmat[5][2][5] = D - ECmat[5][3][4] = E - ECmat[5][3][5] = C124 - ECmat[5][4][0] = F - ECmat[5][4][1] = H - ECmat[5][4][2] = C134 - ECmat[5][4][3] = E - ECmat[5][5][0] = B - ECmat[5][5][1] = C - ECmat[5][5][2] = D - ECmat[5][5][3] = C124 - - - elif(168 <= self.SGN and self.SGN <= 176): # Hexagonal II - LC = 'HII' - ECs= 12 - - A = C111-C222+C112 - B = float(3/4)*C222-float(1/2)*C111-float(1/4)*C112 - C = float(1/2)*C111-float(1/4)*C222-float(1/4)*C112 - D = float(1/2)*(C113-C123) - E = float(1/2)*(C155-C144) - - ECmat[0][0][0] = C111 - ECmat[0][0][1] = C112 - ECmat[0][0][2] = C113 - ECmat[0][0][5] = C116 - ECmat[0][1][0] = C112 - ECmat[0][1][1] = A - ECmat[0][1][2] = C123 - ECmat[0][1][5] = -C116 - ECmat[0][2][0] = C113 - ECmat[0][2][1] = C123 - ECmat[0][2][2] = C133 - ECmat[0][3][3] = C144 - ECmat[0][3][4] = C145 - ECmat[0][4][3] = C145 - ECmat[0][4][4] = C155 - ECmat[0][5][0] = C116 - ECmat[0][5][1] = -C116 - ECmat[0][5][5] = B - - ECmat[1][0][0] = C112 - ECmat[1][0][1] = A - ECmat[1][0][2] = C123 - ECmat[1][1][0] = A - ECmat[1][1][1] = A222 - ECmat[1][1][2] = C113 - ECmat[1][1][5] = C116 - ECmat[1][2][0] = C123 - ECmat[1][2][1] = C113 - ECmat[1][2][2] = C133 - ECmat[1][3][3] = C155 - ECmat[1][3][4] = -C145 - ECmat[1][4][3] = -C145 - ECmat[1][4][4] = C144 - ECmat[1][5][1] = C116 - ECmat[1][5][5] = C - - ECmat[2][0][0] = C113 - ECmat[2][0][1] = C123 - ECmat[2][0][2] = C133 - ECmat[2][1][0] = C123 - ECmat[2][1][1] = C113 - ECmat[2][1][2] = C133 - ECmat[2][2][0] = C133 - ECmat[2][2][1] = C133 - ECmat[2][2][2] = C133 - ECmat[2][3][3] = C344 - ECmat[2][4][4] = C344 - ECmat[2][5][5] = D - - ECmat[3][0][3] = C144 - ECmat[3][1][3] = C155 - ECmat[3][2][3] = C344 - ECmat[3][3][0] = C144 - ECmat[3][3][1] = C155 - ECmat[3][3][2] = C344 - ECmat[3][3][5] = C145 - ECmat[3][4][5] = E - ECmat[3][5][3] = C145 - ECmat[3][5][4] = E - - ECmat[4][0][3] = C155 - ECmat[4][1][3] = C144 - ECmat[4][2][3] = C344 - ECmat[4][3][5] = E - ECmat[4][4][0] = C155 - ECmat[4][4][1] = C144 - ECmat[4][4][2] = C344 - ECmat[4][4][5] = -C145 - ECmat[4][5][3] = E - ECmat[4][5][4] = -C145 - - ECmat[5][0][5] = B - ECmat[5][1][5] = C - ECmat[5][2][5] = D - ECmat[5][3][4] = E - ECmat[5][4][3] = E - ECmat[5][5][0] = B - ECmat[5][5][1] = C - ECmat[5][5][2] = D - ECmat[5][5][5] = -C116 - - elif(177 <= self.SGN and self.SGN <= 194): # Hexagonal I - LC = 'HI' - ECs= 10 - - A = C111-C222+C112 - B = float(3/4)*C222-float(1/2)*C111-float(1/4)*C112 - C = float(1/2)*C111-float(1/4)*C222-float(1/4)*C112 - D = float(1/2)*(C113-C123) - E = float(1/2)*(C155-C144) - - ECmat[0][0][0] = C111 - ECmat[0][0][1] = C112 - ECmat[0][0][2] = C113 - ECmat[0][1][0] = C112 - ECmat[0][1][1] = A - ECmat[0][1][2] = C123 - ECmat[0][2][0] = C113 - ECmat[0][2][1] = C123 - ECmat[0][2][2] = C133 - ECmat[0][3][3] = C144 - ECmat[0][4][4] = C155 - ECmat[0][5][5] = B - - ECmat[1][0][0] = C112 - ECmat[1][0][1] = A - ECmat[1][0][2] = C123 - ECmat[1][1][0] = A - ECmat[1][1][1] = C222 - ECmat[1][1][2] = C113 - ECmat[1][2][0] = C123 - ECmat[1][2][1] = C113 - ECmat[1][2][2] = C133 - ECmat[1][3][3] = C155 - ECmat[1][4][4] = C144 - ECmat[1][5][5] = C - - ECmat[2][0][0] = C113 - ECmat[2][0][1] = C123 - ECmat[2][0][2] = C133 - ECmat[2][1][0] = C123 - ECmat[2][1][1] = C113 - ECmat[2][1][2] = C133 - ECmat[2][2][0] = C133 - ECmat[2][2][1] = C133 - ECmat[2][2][2] = C133 - ECmat[2][3][3] = C344 - ECmat[2][4][4] = C344 - ECmat[2][5][5] = D - - ECmat[3][0][3] = C144 - ECmat[3][1][3] = C155 - ECmat[3][2][3] = C344 - ECmat[3][3][0] = C144 - ECmat[3][3][1] = C155 - ECmat[3][3][2] = C344 - ECmat[3][4][5] = E - ECmat[3][5][4] = E - - ECmat[4][0][4] = C155 - ECmat[4][1][4] = C144 - ECmat[4][2][4] = C344 - ECmat[4][3][5] = E - ECmat[4][4][0] = C155 - ECmat[4][4][1] = C144 - ECmat[4][4][2] = C344 - ECmat[4][5][3] = E - - ECmat[5][0][5] = B - ECmat[5][1][5] = C - ECmat[5][2][5] = D - ECmat[5][3][4] = E - ECmat[5][4][3] = E - ECmat[5][5][0] = B - ECmat[5][5][1] = C - ECmat[5][5][2] = D - - elif(195 <= self.SGN and self.SGN <= 206): # Cubic II - LC = 'CII' - ECs= 8 - - ECmat[0][0][0] = C111 - ECmat[0][0][1] = C112 - ECmat[0][0][2] = C113 - ECmat[0][1][0] = C112 - ECmat[0][1][1] = C113 - ECmat[0][1][2] = C123 - ECmat[0][2][0] = C113 - ECmat[0][2][1] = C123 - ECmat[0][2][2] = C112 - ECmat[0][3][3] = C144 - ECmat[0][4][4] = C155 - ECmat[0][5][5] = C166 - - ECmat[1][0][0] = C112 - ECmat[1][0][1] = C112 - ECmat[1][0][2] = C113 - ECmat[1][1][0] = C112 - ECmat[1][1][1] = C111 - ECmat[1][1][2] = C112 - ECmat[1][2][0] = C123 - ECmat[1][2][1] = C112 - ECmat[1][2][2] = C113 - ECmat[1][3][3] = C166 - ECmat[1][4][4] = C144 - ECmat[1][5][5] = C155 - - ECmat[2][0][0] = C112 - ECmat[2][0][1] = C123 - ECmat[2][0][2] = C112 - ECmat[2][1][0] = C123 - ECmat[2][1][1] = C112 - ECmat[2][1][2] = C112 - ECmat[2][2][0] = C112 - ECmat[2][2][1] = C112 - ECmat[2][2][2] = C111 - ECmat[2][3][3] = C155 - ECmat[2][4][4] = C166 - ECmat[2][5][5] = C144 - - ECmat[3][0][3] = C144 - ECmat[3][1][3] = C155 - ECmat[3][2][3] = C155 - ECmat[3][3][0] = C144 - ECmat[3][3][1] = C155 - ECmat[3][3][2] = C155 - ECmat[3][4][5] = C456 - ECmat[3][5][4] = C456 - - ECmat[4][0][3] = C155 - ECmat[4][1][3] = C144 - ECmat[4][2][3] = C155 - ECmat[4][3][5] = C456 - ECmat[4][4][0] = C155 - ECmat[4][4][1] = C144 - ECmat[4][4][2] = C155 - ECmat[4][5][3] = C456 - - ECmat[5][0][5] = C155 - ECmat[5][1][5] = C155 - ECmat[5][2][5] = C144 - ECmat[5][3][4] = C456 - ECmat[5][4][3] = C456 - ECmat[5][5][0] = C155 - ECmat[5][5][1] = C155 - ECmat[5][5][2] = C144 - - elif(207 <= self.SGN and self.SGN <= 230): # Cubic I - LC = 'CI' - ECs= 6 - - ECmat[0][0][0] = C111 - ECmat[0][0][1] = C112 - ECmat[0][0][2] = C112 - ECmat[0][1][0] = C112 - ECmat[0][1][1] = C112 - ECmat[0][1][2] = C113 - ECmat[0][2][0] = C112 - ECmat[0][2][1] = C123 - ECmat[0][2][2] = C112 - ECmat[0][3][3] = C144 - ECmat[0][4][4] = C155 - ECmat[0][5][5] = C155 - - ECmat[1][0][0] = C112 - ECmat[1][0][1] = C112 - ECmat[1][0][2] = C123 - ECmat[1][1][0] = C112 - ECmat[1][1][1] = C111 - ECmat[1][1][2] = C112 - ECmat[1][2][0] = C123 - ECmat[1][2][1] = C112 - ECmat[1][2][2] = C112 - ECmat[1][3][3] = C155 - ECmat[1][4][4] = C144 - ECmat[1][5][5] = C155 - - ECmat[2][0][0] = C112 - ECmat[2][0][1] = C123 - ECmat[2][0][2] = C112 - ECmat[2][1][0] = C123 - ECmat[2][1][1] = C112 - ECmat[2][1][2] = C112 - ECmat[2][2][0] = C112 - ECmat[2][2][1] = C112 - ECmat[2][2][2] = C111 - ECmat[2][3][3] = C155 - ECmat[2][4][4] = C155 - ECmat[2][5][5] = C144 - - ECmat[3][0][3] = C144 - ECmat[3][1][3] = C155 - ECmat[3][2][3] = C155 - ECmat[3][3][0] = C144 - ECmat[3][3][1] = C155 - ECmat[3][3][2] = C155 - ECmat[3][4][5] = C456 - ECmat[3][5][4] = C456 - - ECmat[4][0][3] = C155 - ECmat[4][1][3] = C144 - ECmat[4][2][3] = C155 - ECmat[4][3][5] = C456 - ECmat[4][4][0] = C155 - ECmat[4][4][1] = C144 - ECmat[4][4][2] = C155 - ECmat[4][5][3] = C456 - - ECmat[5][0][5] = C155 - ECmat[5][1][5] = C155 - ECmat[5][2][5] = C144 - ECmat[5][3][4] = C456 - ECmat[5][4][3] = C456 - ECmat[5][5][0] = C155 - ECmat[5][5][1] = C155 - ECmat[5][5][2] = C144 - - elasticSIndex = backend.openSection("x_elastic_section_strain_diagrams") - backend.addValue("x_elastic_strain_diagram_type", "energy") - backend.addValue("x_elastic_strain_diagram_number_of_eta", len(eta)) - backend.addValue("x_elastic_strain_diagram_eta_values", eta) - backend.addValue("x_elastic_strain_diagram_values", energy) - backend.closeSection("x_elastic_section_strain_diagrams", elasticSIndex) - backend.addValue('x_elastic_3rd_order_constants_matrix',ECmat) - backend.closeSection("section_single_configuration_calculation", elasticGIndex) - backend.addValue("x_elastic_deformation_types", defTyp) - backend.addValue("x_elastic_number_of_deformations", defNum) - if meth == "Energy": - elasticPIndex = backend.openSection("x_elastic_section_fitting_parameters") - backend.addValue("x_elastic_fitting_parameters_eta", self.etaEC) - backend.addValue("x_elastic_fitting_parameters_polynomial_order", self.fitEC) - backend.closeSection("x_elastic_section_fitting_parameters", elasticPIndex) - else: - pass - - def onClose_section_single_configuration_calculation(self, backend, gIndex, section): -# logging.error("BASE onClose_section_single_configuration_calculation") - backend.addValue('single_configuration_to_calculation_method_ref', self.secMethodIndex) - backend.addValue('single_configuration_calculation_to_system_ref', self.secSystemIndex) - -mainFileDescription = \ - SM(name = 'root', - weak = False, - startReStr = "", - subMatchers = [ - SM(name = 'input', - startReStr = r"\s*Order of elastic constants\s*=\s*(?P<x_elastic_elastic_constant_order>[0-9]+)", - repeats = False, - required = False, - forwardMatch = False, - sections = ['section_run', 'section_method'], - subMatchers = [ - SM(r"\s*Method of calculation\s*=\s*(?P<x_elastic_calculation_method>[-a-zA-Z]+)"), - SM(r"\s*DFT code name\s*=\s*(?P<x_elastic_code>[-a-zA-Z]+)"), - SM(name = 'system', - startReStr = r"\s*Space-group number\s*=\s*(?P<x_elastic_space_group_number>[0-9]+)", - sections = ['section_system'], - subMatchers = [ - SM(r"\s*Volume of equilibrium unit cell\s*=\s*(?P<x_elastic_unit_cell_volume__bohr3>[-0-9.]+)\s*\[a.u\^3\]") - ]), - SM(r"\s*Maximum Lagrangian strain\s*=\s*(?P<x_elastic_max_lagrangian_strain>[0-9.]+)"), - SM(r"\s*Number of distorted structures\s*=\s*(?P<x_elastic_number_of_distorted_structures>[0-9]+)") - ] ) - ]) - - -parserInfo = { - "name": "elastic_parser", - "version": "1.0" -} - - -class ElasticParser(): - """ A proper class envolop for running this parser from within python. """ - def __init__(self, backend, **kwargs): - self.backend_factory = backend - - def parse(self, mainfile): - from unittest.mock import patch - logging.info('elastic parser started') - logging.getLogger('nomadcore').setLevel(logging.WARNING) - backend = self.backend_factory("elastic.nomadmetainfo.json") - with patch.object(sys, 'argv', ['<exe>', '--uri', 'nmd://uri', mainfile]): - mainFunction( - mainFileDescription, - None, - parserInfo, - superContext=SampleContext(), - superBackend=backend) - - return backend -- GitLab