diff --git a/elasticparser/elastic_parser.py b/elasticparser/elastic_parser.py
index 67311ba54122d1a952d9e80d5e7b1f0b8956182a..ac236ad77a16cc08470dd65c61a9d9c79a333b9c 100644
--- a/elasticparser/elastic_parser.py
+++ b/elasticparser/elastic_parser.py
@@ -2,7 +2,6 @@ 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
@@ -12,10 +11,6 @@ from elasticparser.metainfo.elastic import x_elastic_section_strain_diagrams,\
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):
@@ -40,7 +35,7 @@ class ElasticParserInterface:
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
+ sec_strain_diagram.x_elastic_strain_diagram_values = energy
poly_fit_2 = int((n_strains - 1) / 2)
@@ -60,8 +55,7 @@ class ElasticParserInterface:
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
+ sec_strain_diagram.x_elastic_strain_diagram_values = energy_fit[diagram_type][1][fit_order]
elif method == 'stress':
strain, stress = self.properties.get_strain_stress()
@@ -92,8 +86,7 @@ class ElasticParserInterface:
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
+ sec_strain_diagram.x_elastic_strain_diagram_values = np.array(stress_fit[diagram_type][si][1][fit_order])
def parse_elastic_constant(self):
sec_scc = self.archive.section_run[-1].section_single_configuration_calculation[-1]
@@ -104,26 +97,26 @@ class ElasticParserInterface:
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_2nd_order_constants_matrix = matrices['elastic_constant']
+ sec_scc.x_elastic_2nd_order_constants_compliance_matrix = matrices['compliance']
- 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_Voigt_bulk_modulus = moduli.get('B_V', moduli.get('K_V'))
+ sec_scc.x_elastic_Voigt_shear_modulus = moduli['G_V']
- 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_Reuss_bulk_modulus = moduli.get('B_R', moduli.get('K_R'))
+ sec_scc.x_elastic_Reuss_shear_modulus = moduli['G_R']
- 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_Hill_bulk_modulus = moduli.get('B_H', moduli.get('K_H'))
+ sec_scc.x_elastic_Hill_shear_modulus = moduli['G_H']
- sec_scc.x_elastic_Voigt_Young_modulus = moduli['E_V'] * giga
+ sec_scc.x_elastic_Voigt_Young_modulus = moduli['E_V']
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_Young_modulus = moduli['E_R']
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_Young_modulus = moduli['E_H']
sec_scc.x_elastic_Hill_Poisson_ratio = moduli['nu_H']
- sec_scc.x_elastic_eigenvalues = eigenvalues * giga
+ sec_scc.x_elastic_eigenvalues = eigenvalues
elif order == 3:
elastic_constant = self.properties.get_elastic_constants_order3()
@@ -141,11 +134,10 @@ class ElasticParserInterface:
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.atom_positions = positions
+ sec_system.simulation_cell = cell
sec_system.configuration_periodic_dimensions = [True, True, True]
sec_system.x_elastic_space_group_number = self.properties.info['space_group_number']
sec_system.x_elastic_unit_cell_volume = volume
diff --git a/elasticparser/elastic_properties.py b/elasticparser/elastic_properties.py
index 4047b063c99e135ea68a5cee3ec2c36c701695ce..aab6dfa5cd9c32f68b9758754af4b4d0126ef818 100644
--- a/elasticparser/elastic_properties.py
+++ b/elasticparser/elastic_properties.py
@@ -1,4 +1,5 @@
import os
+import pint
import numpy as np
from ase import Atoms
@@ -32,7 +33,7 @@ class ElasticProperties:
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*'),
+ Quantity('equilibrium_volume', r'\s*Volume of equilibrium unit cell\s*=\s*([0-9.]+)\s*', unit='angstrom ** 3'),
Quantity('max_strain', r'\s*Maximum Lagrangian strain\s*=\s*([0-9.]+)'),
Quantity('n_strains', r'\s*Number of distorted structures\s*=\s*([0-9]+)')
]
@@ -107,7 +108,12 @@ class ElasticProperties:
structure = Atoms(cell=cellpar, scaled_positions=pos, symbols=sym, pbc=True)
- return sym, structure.get_positions(), list(structure.get_cell())
+ positions = structure.get_positions()
+ positions = pint.Quantity(positions, 'angstrom')
+ cell = structure.get_cell()
+ cell = pint.Quantity(cell, 'angstrom')
+
+ return sym, positions, cell
def get_strain_energy(self):
strains, energies = [], []
@@ -122,6 +128,11 @@ class ElasticProperties:
strains.append(data[0])
energies.append(data[1])
+ energies = pint.Quantity(energies, 'hartree')
+ # the peculiarity of the x_elastic_strain_diagram_values metainfo that it does
+ # not have the energy unit
+ energies = energies.to('J').magnitude
+
return strains, energies
def get_strain_stress(self):
@@ -206,7 +217,7 @@ class ElasticProperties:
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
+ return [eta, val]
def get_energy_fit(self):
energy_fit = dict()
@@ -216,10 +227,14 @@ class ElasticProperties:
if result is None:
continue
+ result[1] = {
+ key: pint.Quantity(val, 'GPa').to('Pa').magnitude for key, val in result[1].items()}
energy_fit['d2e'] = result
result = self._get_fit('Energy-vs-Strain', 'CVe.dat')
if result is not None:
+ result[1] = {
+ key: pint.Quantity(val, 'hartree').to('J').magnitude for key, val in result[1].items()}
energy_fit['cross-validation'] = result
return energy_fit
@@ -232,10 +247,12 @@ class ElasticProperties:
for strain_index in range(1, 7):
result = self._get_fit('Stress-vs-Strain', '%d_dS.dat' % strain_index)
if result is not None:
+ result[1] = {key: pint.Quantity(val, 'GPa') for key, val in result[1].items()}
stress_fit['dtn'][strain_index - 1] = result
result = self._get_fit('Stress-vs-Strain', '%d_CVe.dat' % strain_index)
if result is not None:
+ result[1] = {key: pint.Quantity(val, 'hartree') for key, val in result[1].items()}
stress_fit['cross-validation'][strain_index - 1] = result
return stress_fit
@@ -301,21 +318,25 @@ class ElasticProperties:
str_operation=reshape, dtype=str),
Quantity(
'elastic_constant', r'Elastic constant[\s\S]+in GPa\s*:\s*\n\n([\-\d\.\s\n]+)\n',
- str_operation=reshape, dtype=float),
+ str_operation=reshape, dtype=float, unit='GPa'),
Quantity(
'compliance', r'Elastic compliance[\s\S]+in 1/GPa\s*:\s*\n\n([\-\d\.\s\n]+)\n',
- str_operation=reshape, dtype=float)]
+ str_operation=reshape, dtype=float, unit='1/GPa')]
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]
+ for name in modulus_names:
+ unit = None
+ if name[0:2] in ('B_', 'K_', 'G_', 'E_'):
+ unit = 'GPa'
+ quantities.append(
+ Quantity(name, r'%s\s*=\s*([-+\d+\.]+)\s*' % name, unit=unit))
quantities += [
Quantity(
- 'eigenvalues', r'Eigenvalues of elastic constant \(stiffness\) matrix:\s*\n+([\-\d\.\n\s]+)\n')]
+ 'eigenvalues', r'Eigenvalues of elastic constant \(stiffness\) matrix:\s*\n+([\-\d\.\n\s]+)\n', unit='GPa')]
parser = UnstructuredTextFileParser(path, quantities)
diff --git a/elasticparser/metainfo/elastic.py b/elasticparser/metainfo/elastic.py
index 1928e0cd305a9aa21ea5c23b8af4dc48834c1486..ba5a444b192e4ef35b40c0578c1fbfe6ca12371f 100644
--- a/elasticparser/metainfo/elastic.py
+++ b/elasticparser/metainfo/elastic.py
@@ -330,6 +330,7 @@ class section_system(public.section_system):
x_elastic_unit_cell_volume = Quantity(
type=np.dtype(np.float64),
shape=[],
+ unit='m ** 3',
description='''
Volume of the equilibrium unit cell
''',