Skip to content
Snippets Groups Projects
Commit c0a152ed authored by Alvin Noe Ladines's avatar Alvin Noe Ladines
Browse files

Switched to use nomad metainfo

parent f41d95aa
Branches
No related tags found
No related merge requests found
...@@ -12,4 +12,21 @@ ...@@ -12,4 +12,21 @@
# See the License for the specific language governing permissions and # See the License for the specific language governing permissions and
# limitations under the License. # 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()
# 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
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']
# 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)
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
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
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment