# -*- coding: utf-8 -*-
from aiida.orm.calculation.job.quantumespresso.cp import CpCalculation
from aiida.parsers.plugins.quantumespresso.basic_raw_parser_cp import (
QEOutputParsingError, parse_cp_traj_stanzas, parse_cp_raw_output)
from aiida.parsers.plugins.quantumespresso.constants import (bohr_to_ang,
timeau_to_sec, hartree_to_ev)
from aiida.orm.data.parameter import ParameterData
from aiida.orm.data.structure import StructureData
from aiida.orm.data.folder import FolderData
from aiida.parsers.parser import Parser
from aiida.common.datastructures import calc_states
from aiida.orm.data.array.trajectory import TrajectoryData
import numpy
__copyright__ = u"Copyright (c), 2015, ECOLE POLYTECHNIQUE FEDERALE DE LAUSANNE (Theory and Simulation of Materials (THEOS) and National Centre for Computational Design and Discovery of Novel Materials (NCCR MARVEL)), Switzerland and ROBERT BOSCH LLC, USA. All rights reserved."
__license__ = "MIT license, see LICENSE.txt file"
__version__ = "0.5.0"
__contributors__ = "Andrea Cepellotti, Giovanni Pizzi, Martin Uhrin, Nicolas Mounet"
[docs]class CpParser(Parser):
"""
This class is the implementation of the Parser class for Cp.
"""
def __init__(self, calc):
"""
Initialize the instance of CpParser
:param calculation: calculation object.
"""
# check for valid input
if not isinstance(calc, CpCalculation):
raise QEOutputParsingError("Input calc must be a CpCalculation")
super(CpParser, self).__init__(calc)
[docs] def parse_with_retrieved(self, retrieved):
"""
Receives in input a dictionary of retrieved nodes.
Does all the logic here.
"""
from aiida.common.exceptions import InvalidOperation
import os, copy
import numpy # TrajectoryData also uses numpy arrays
successful = True
# check if I'm not to overwrite anything
state = self._calc.get_state()
if state != calc_states.PARSING:
raise InvalidOperation("Calculation not in {} state"
.format(calc_states.PARSING))
# get the input structure
input_structure = self._calc.inp.structure
# load the input dictionary
# TODO: pass this input_dict to the parser. It might need it.
input_dict = self._calc.inp.parameters.get_dict()
# Check that the retrieved folder is there
try:
out_folder = retrieved[self._calc._get_linkname_retrieved()]
except KeyError:
self.logger.error("No retrieved folder found")
return False, ()
# check what is inside the folder
list_of_files = out_folder.get_folder_list()
# at least the stdout should exist
if not self._calc._OUTPUT_FILE_NAME in list_of_files:
successful = False
new_nodes_tuple = ()
self.logger.error("Standard output not found")
return successful, new_nodes_tuple
# if there is something more, I note it down, so to call the raw parser
# with the right options
# look for xml
out_file = out_folder.get_abs_path(self._calc._OUTPUT_FILE_NAME)
xml_file = None
if self._calc._DATAFILE_XML_BASENAME in list_of_files:
xml_file = out_folder.get_abs_path(self._calc._DATAFILE_XML_BASENAME)
xml_counter_file = None
if self._calc._FILE_XML_PRINT_COUNTER in list_of_files:
xml_counter_file = out_folder.get_abs_path(
self._calc._FILE_XML_PRINT_COUNTER)
parsing_args = [out_file, xml_file, xml_counter_file]
# call the raw parsing function
out_dict, raw_successful = parse_cp_raw_output(*parsing_args)
successful = True if raw_successful else False
# parse the trajectory. Units in Angstrom, picoseconds and eV.
# append everthing in the temporary dictionary raw_trajectory
expected_configs = None
raw_trajectory = {}
evp_keys = ['electronic_kinetic_energy', 'cell_temperature', 'ionic_temperature',
'scf_total_energy', 'enthalpy', 'enthalpy_plus_kinetic',
'energy_constant_motion', 'volume', 'pressure']
pos_vel_keys = ['cells', 'positions', 'times', 'velocities']
# set a default null values
# Now prepare the reordering, as filex in the xml are ordered
reordering = self._generate_sites_ordering(out_dict['species'],
out_dict['atoms'])
# =============== POSITIONS trajectory ============================
try:
with open(out_folder.get_abs_path(
'{}.pos'.format(self._calc._PREFIX))) as posfile:
pos_data = [l.split() for l in posfile]
# POSITIONS stored in angstrom
traj_data = parse_cp_traj_stanzas(num_elements=out_dict['number_of_atoms'],
splitlines=pos_data,
prepend_name='positions_traj',
rescale=bohr_to_ang)
# here initialize the dictionary. If the parsing of positions fails, though, I don't have anything
# out of the CP dynamics. Therefore, the calculation status is set to FAILED.
raw_trajectory['positions_ordered'] = self._get_reordered_array(traj_data['positions_traj_data'],
reordering)
raw_trajectory['times'] = numpy.array(traj_data['positions_traj_times'])
except IOError:
out_dict['warnings'].append("Unable to open the POS file... skipping.")
successful = False
except Exception as e:
out_dict['warnings'].append("Error parsing POS file ({}). Skipping file."
.format(e.message))
successful = False
# =============== CELL trajectory ============================
try:
with open(os.path.join(out_folder.get_abs_path('.'),
'{}.cel'.format(self._calc._PREFIX))) as celfile:
cel_data = [l.split() for l in celfile]
traj_data = parse_cp_traj_stanzas(num_elements=3,
splitlines=cel_data,
prepend_name='cell_traj',
rescale=bohr_to_ang)
raw_trajectory['cells'] = numpy.array(traj_data['cell_traj_data'])
except IOError:
out_dict['warnings'].append("Unable to open the CEL file... skipping.")
except Exception as e:
out_dict['warnings'].append("Error parsing CEL file ({}). Skipping file."
.format(e.message))
# =============== VELOCITIES trajectory ============================
try:
with open(os.path.join(out_folder.get_abs_path('.'),
'{}.vel'.format(self._calc._PREFIX))) as velfile:
vel_data = [l.split() for l in velfile]
traj_data = parse_cp_traj_stanzas(num_elements=out_dict['number_of_atoms'],
splitlines=vel_data,
prepend_name='velocities_traj',
rescale=bohr_to_ang / timeau_to_sec * 10 ** 12) # velocities in ang/ps,
raw_trajectory['velocities_ordered'] = self._get_reordered_array(traj_data['velocities_traj_data'],
reordering)
except IOError:
out_dict['warnings'].append("Unable to open the VEL file... skipping.")
except Exception as e:
out_dict['warnings'].append("Error parsing VEL file ({}). Skipping file."
.format(e.message))
# =============== EVP trajectory ============================
try:
matrix = numpy.genfromtxt(os.path.join(out_folder.get_abs_path('.'),
'{}.evp'.format(self._calc._PREFIX)))
# there might be a different format if the matrix has one row only
try:
matrix.shape[1]
except IndexError:
matrix = numpy.array(numpy.matrix(matrix))
raw_trajectory['steps'] = numpy.array(matrix[:, 0], dtype=int)
raw_trajectory['electronic_kinetic_energy'] = matrix[:, 1] * hartree_to_ev # EKINC, eV
raw_trajectory['cell_temperature'] = matrix[:, 2] # TEMPH, K
raw_trajectory['ionic_temperature'] = matrix[:, 3] # TEMPP, K
raw_trajectory['scf_total_energy'] = matrix[:, 4] * hartree_to_ev # ETOT, eV
raw_trajectory['enthalpy'] = matrix[:, 5] * hartree_to_ev # ENTHAL, eV
raw_trajectory['enthalpy_plus_kinetic'] = matrix[:, 6] * hartree_to_ev # ECONS, eV
raw_trajectory['energy_constant_motion'] = matrix[:, 7] * hartree_to_ev # ECONT, eV
raw_trajectory['volume'] = matrix[:, 8] * (bohr_to_ang ** 3) # volume, angstrom^3
raw_trajectory['pressure'] = matrix[:, 9] # out_press, GPa
except Exception as e:
out_dict['warnings'].append("Error parsing EVP file ({}). Skipping file.".format(e.message))
except IOError:
out_dict['warnings'].append("Unable to open the EVP file... skipping.")
# get the symbols from the input
# TODO: I should have kinds in TrajectoryData
raw_trajectory['symbols'] = numpy.array([str(i.kind_name) for i in input_structure.sites])
traj = TrajectoryData()
traj.set_trajectory(steps=raw_trajectory['steps'],
cells=raw_trajectory['cells'],
symbols=raw_trajectory['symbols'],
positions=raw_trajectory['positions_ordered'],
times=raw_trajectory['times'],
velocities=raw_trajectory['velocities_ordered'],
)
for this_name in evp_keys:
traj.set_array(this_name, raw_trajectory[this_name])
new_nodes_list = [(self.get_linkname_trajectory(), traj)]
# convert the dictionary into an AiiDA object
output_params = ParameterData(dict=out_dict)
# save it into db
new_nodes_list.append((self.get_linkname_outparams(), output_params))
return successful, new_nodes_list
[docs] def get_linkname_trajectory(self):
"""
Returns the name of the link to the output_structure (None if not present)
"""
return 'output_trajectory'
def _generate_sites_ordering(self, raw_species, raw_atoms):
"""
take the positions of xml and from file.pos of the LAST step and compare them
"""
# Examples in the comments are for species [Ba, O, Ti]
# and atoms [Ba, Ti, O, O, O]
# Dictionary to associate the species name to the idx
# Example: {'Ba': 1, 'O': 2, 'Ti': 3}
species_dict = {name: idx for idx, name in zip(raw_species['index'],
raw_species['type'])}
# List of the indices of the specie associated to each atom,
# in the order specified in input
# Example: (1,3,2,2,2)
atoms_species_idx = [species_dict[a[0]] for a in raw_atoms]
# I also attach the current position; important to convert to a list
# Otherwise the iterator can be looped on only once!
# Example: ((0,1),(1,3),(2,2),(3,2),(4,2))
ref_atom_list = list(enumerate(atoms_species_idx))
new_order_tmp = []
# I reorder the atoms, first by specie, then in their order
# This is the order used in output by CP!!
# Example: ((0,1),(2,2),(3,2),(4,2),(1,3))
for specie_idx in sorted(raw_species['index']):
for elem in ref_atom_list:
if elem[1] == specie_idx:
new_order_tmp.append(elem)
# This is the new order that is printed in CP:
# e.g. reordering[2] is the index of the atom, in the input
# list of atoms, that is printed in position 2 (0-based, so the
# third atom) in the CP output files.
# Example: [0,2,3,4,1]
reordering = [_[0] for _ in new_order_tmp]
# I now need the inverse reordering, to put back in place
# from the output ordering to the input one!
# Example: [0,4,1,2,3]
# Because in the final list (Ba, O, O, O, Ti)
# the first atom Ba in the input is atom 0 in the CP output (the first),
# the second atom Ti in the input is atom 4 (the fifth) in the CP output,
# and so on
sorted_indexed_reordering = sorted([(_[1], _[0]) for _ in
enumerate(reordering)])
reordering_inverse = [_[1] for _ in sorted_indexed_reordering]
return reordering_inverse
def _get_reordered_list(self, origlist, reordering):
"""
Given a list to reorder, a list of integer positions with the new
order, return the reordered list.
"""
return [origlist[e] for e in reordering]
def _get_reordered_array(self, input, reordering):
return numpy.array([self._get_reordered_list(i, reordering) for i in input])