Commit 42a6cb6f authored by Mikkel Strange's avatar Mikkel Strange

forces and energies extracted from Analysis in ATK

parent dcc7da77
......@@ -7,25 +7,75 @@ import re
class Reader:
def __init__(self, fname):
self.atoms_x = {} # like {'gID0001': Atoms('H2')}
self.calculator_x = {} # like {'gID0001': LCAOCalculator}
self.atoms_x = {} # like {'gID0001': Atoms('H2')}
self.calculator_x = {} # like {'gID0001': LCAOCalculator}
self.conf_names = None
self.calc_names = None
self.f = netcdf_file(fname, 'r', mmap=True)
self.initialize()
gids = self.calc_names.keys()
for gid in gids:
self.atk_version = self.f.version[:].decode('utf-8').split()[-1]
self.read_names()
for gid in self.calc_names.keys():
conf_name = self.conf_names[gid]
calc_name = self.calc_names[gid]
fpt = self.finger_print_table[gid] # req for for analysis parsing
self.atoms_x[gid] = parse_configuration(self.f, conf_name)
self.calculator_x[gid] = parse_calculator(self.f, calc_name)
def initialize(self):
# setup hamilton, results, wave_functions attr
# to calculator to hold results (as in GPAW)
p_etot_fp = r"^(?P<pref>TotalEnergy_gID[0-9]+)_finger_print"
p_forc_fp = r"^(?P<pref>Forces_gID[0-9]+)_finger_print"
for gid in self.calc_names.keys(): # loop over calculators
calc = self.calculator_x[gid]
calc._hamiltonian = type('Hamiltonian', (object,), {})()
calc._results = type('Results', (object,), {})()
calc._wave_functions = type('WaveFunctions', (object,), {})()
fp = self.finger_print_table[gid]
v = self.f.variables
h = calc._hamiltonian
r = calc._results
wf = calc._wave_functions
for name in self.f.variables.keys():
# calc.hamiltonian.e_x, x=e_tot, e_kin, ...
# TotalEnergy object in ATK
m_fp_etot = re.search(p_etot_fp, name)
if m_fp_etot is not None:
pref = m_fp_etot.group('pref')
fp2 = v[pref + '_finger_print'].data[:].copy()
fp2 = fp2.tostring().decode('utf-8')
pref += '_component_'
if fp == fp2:
h.e_kin = v[pref + 'Kinetic'].data[:].copy()[0]
h.e_coulomb = v[pref +
'Electrostatic'].data[:].copy()[0]
h.e_xc = v[pref +
'Exchange-Correlation'].data[:].copy()[0]
h.e_entropy = v[pref +
'Entropy-Term'].data[:].copy()[0]
h.e_external = v[pref +
'External-Field'].data[:].copy()[0]
h.e_total_free = (h.e_kin + h.e_coulomb + h.e_xc +
h.e_entropy + h.e_external)
continue
# calc.results.forces
m_fp_forc = re.search(p_forc_fp, name)
if m_fp_forc is not None:
pref = m_fp_forc.group('pref')
fp2 = v[pref + '_finger_print'].data[:].copy()
fp2 = fp2.tostring().decode('utf-8')
if fp == fp2:
r.forces = v[pref +
'_atom_resolved_forces'].data[:].copy()
continue
def read_names(self):
"""Read the names of the variables in the netcdf file for
configurations and calculators and setup
the finger print table which maps between calculated
quantities and configurations.
"""
self.atk_version = self.f.version[:].decode('utf-8').split()[-1]
self._names = self.f._names[:].decode('utf-8').split(';')
self.conf_names = self._read_configuration_names()
self.calc_names = self._read_calculator_names()
self.finger_print_table = self._read_finger_print_table()
......@@ -88,5 +138,5 @@ if __name__ == '__main__':
import sys
r = Reader(sys.argv[1])
for key, value in r.atoms_x.items():
print(key,value)
print(key, value)
print(r.get_atoms(0))
......@@ -46,7 +46,7 @@ def parse(filename):
r = Reader(filename) # Reader(filename)
index = 0 # need to loop over index at some point if more that one conf per
# file
r.calculator = r.get_calculator(index)
r.c = r.get_calculator(index)
r.atoms = r.get_atoms(index)
p.startedParsingSession(filename, parser_info)
with o(p, 'section_run'):
......@@ -81,12 +81,12 @@ def parse(filename):
# c(r.convergence.scf_energy, 'eV')) # eV / electron
p.addValue('smearing_kind', 'fermi')
p.addRealValue('smearing_width',
c(r.calculator.numerical_accuracy_parameters.\
c(r.c.numerical_accuracy_parameters.\
electron_temperature, 'K'))
p.addRealValue('total_charge', r.calculator.charge)
p.addRealValue('total_charge', r.c.charge)
with o(p, 'section_XC_functionals'):
p.addValue('XC_functional_name',
get_libxc_name(r.calculator.exchange_correlation))
get_libxc_name(r.c.exchange_correlation))
with o(p, 'section_single_configuration_calculation'):
p.addValue('single_configuration_calculation_to_system_ref',
system_gid)
......@@ -97,31 +97,31 @@ def parse(filename):
# p.addRealValue('energy_total',
# c(r.hamiltonian.e_tot_extrapolated, 'eV'))
p.addRealValue('energy_free',
c(r.hamiltonian.e_tot, 'eV'))
p.addRealValue('energy_XC', c(r.hamiltonian.e_xc, 'eV'))
c(r.c._hamiltonian.e_total_free, 'eV'))
p.addRealValue('energy_XC', c(r.c._hamiltonian.e_xc, 'eV'))
p.addRealValue('electronic_kinetic_energy',
c(r.hamiltonian.e_kin, 'eV'))
c(r.c._hamiltonian.e_kin, 'eV'))
p.addRealValue('energy_correction_entropy',
c(r.hamiltonian.e_S, 'eV'))
c(r.c._hamiltonian.e_entropy, 'eV'))
# p.addRealValue('energy_reference_fermi',
# c(r.occupations.fermilevel, 'eV'))
if hasattr(r.results, 'forces'):
if hasattr(r.c._results, 'forces'):
p.addArrayValues('atom_forces_free_raw',
c(r.results.forces, 'eV/angstrom'))
c(r.c._results.forces, 'eV/angstrom'))
#if hasattr(r.results, 'magmoms'):
# p.addArrayValues('x_gpaw_magnetic_moments',
# r.results.magmoms)
# p.addRealValue('x_atk_spin_Sz', r.results.magmoms.sum() / 2.0)
if hasattr(r.wave_functions, 'eigenvalues'):
if hasattr(r.c._wave_functions, 'eigenvalues'):
with o(p, 'section_eigenvalues'):
p.addValue('eigenvalues_kind', 'normal')
p.addArrayValues('eigenvalues_values',
c(r.wave_functions.eigenvalues, 'eV'))
c(r.c._wave_functions.eigenvalues, 'eV'))
p.addArrayValues('eigenvalues_occupation',
r.wave_functions.occupations)
r.c._wave_functions.occupations)
p.addArrayValues('eigenvalues_kpoints',
r.wave_functions.ibz_kpts)
if hasattr(r.wave_functions, 'band_paths'):
r.c._wave_functions.ibz_kpts)
if hasattr(r.c._wave_functions, 'band_paths'):
with o(p, 'section_k_band'):
for band_path in r.wave_functions.band_paths:
with o(p, 'section_k_band_segment'):
......
......@@ -61,7 +61,7 @@ def parse_calculator(fd, calcname):
"""
code = fd.variables[calcname].data[:].copy()
code = code.tostring().decode("utf-8")
if 1:
if 0:
print(code)
#s = re.search('\s*(?P<name>[0-9a-zA-Z_]+)\s*=\s*LCAOCalculator\(', code)
#name = s.group('name')
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment