Added support for phonony db, switched to metainfo

phonopyparser/FHIaims.py

+# FHIaims.py - IO routines for phonopy-FHI-aims
+# methods compatible with the corresponding ones from ase.io.aims
+# only minimal subset of functionality required within phonopy context is implemented
+#
+# Copyright (C) 2009-2011 Joerg Meyer (jm)
+# All rights reserved.
+#
+# This file is part of phonopy.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions
+# are met:
+#
+# * Redistributions of source code must retain the above copyright
+#   notice, this list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright
+#   notice, this list of conditions and the following disclaimer in
+#   the documentation and/or other materials provided with the
+#   distribution.
+#
+# * Neither the name of the phonopy project nor the names of its
+#   contributors may be used to endorse or promote products derived
+#   from this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+# COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+
+from phonopy.structure.atoms import PhonopyAtoms
+
+def read_aims(filename):
+    """Method to read FHI-aims geometry files in phonopy context."""
+
+    lines = open(filename, 'r').readlines()
+
+    cell = []
+    is_frac = []
+    positions = []
+    symbols = []
+    magmoms = []
+    for line in lines:
+        fields = line.split()
+        if not len(fields):
+            continue
+        if fields[0] == "lattice_vector":
+            vec = [float(x) for x in fields[1:4]]
+            cell.append(vec)
+        elif fields[0][0:4] == "atom":
+            if fields[0] == "atom":
+                frac = False
+            elif fields[0] == "atom_frac":
+                frac = True
+            pos = [float(x) for x in fields[1:4]]
+            sym = fields[4]
+            is_frac.append(frac)
+            positions.append(pos)
+            symbols.append(sym)
+            magmoms.append(None)
+        # implicitly assuming that initial_moments line adhere to FHI-aims geometry.in specification,
+        # i.e. two subsequent initial_moments lines do not occur
+        # if they do, the value specified in the last line is taken here - without any warning
+        elif fields[0] == "initial_moment":
+            magmoms[-1] = float(fields[1])
+
+    for (n,frac) in enumerate(is_frac):
+        if frac:
+            pos = [ sum( [ positions[n][l] * cell[l][i] for l in range(3) ] ) for i in range(3) ]
+            positions[n] = pos
+    if None in magmoms:
+        atoms = PhonopyAtoms(cell=cell, symbols=symbols, positions=positions)
+    else:
+        atoms = PhonopyAtoms(cell=cell, symbols=symbols, positions=positions, magmoms=magmoms)
+
+    return atoms
+
+
+def write_aims(filename, atoms):
+    """Method to write FHI-aims geometry files in phonopy context."""
+
+    lines = ""
+    lines += "# geometry.in for FHI-aims \n"
+    lines += "# | generated by phonopy.FHIaims.write_aims() \n"
+
+    lattice_vector_line = "lattice_vector " + "%16.16f "*3 + "\n"
+    for vec in atoms.get_cell():
+        lines += lattice_vector_line % tuple(vec)
+
+    N = atoms.get_number_of_atoms()
+
+    atom_line = "atom " + "%16.16f "*3 + "%s \n"
+    positions = atoms.get_positions()
+    symbols = atoms.get_chemical_symbols()
+
+    initial_moment_line = "initial_moment %16.6f\n"
+    magmoms = atoms.get_magnetic_moments()
+
+    for n in range(N):
+        lines += atom_line % (tuple(positions[n]) + (symbols[n],))
+        if magmoms is not None:
+            lines += initial_moment_line % magmoms[n]
+
+    f = open(filename, 'w')
+    f.write(lines)
+    f.close()
+
+
+
+class Atoms_with_forces(PhonopyAtoms):
+    """ Hack to phonopy.atoms to maintain ASE compatibility also for forces."""
+
+    def get_forces(self):
+        return self.forces
+
+
+def read_aims_output(filename):
+    """ Read FHI-aims output and 
+        return geometry, energy and forces from last self-consistency iteration"""
+
+    lines = open(filename, 'r').readlines()
+
+    l = 0
+    N = 0
+    while l < len(lines):
+        line = lines[l]
+        if "Number of atoms" in line:
+            N = int(line.split()[5])
+        elif "| Unit cell:" in line:
+            cell = []
+            for i in range(3):
+                l += 1
+                vec = [float(x) for x in lines[l].split()[1:4]]
+                cell.append(vec)
+        elif ("Atomic structure:" in line) or ("Updated atomic structure:" in line):
+            if "Atomic structure:" in line:
+                i_sym = 3
+                i_pos_min = 4 ; i_pos_max = 7
+            elif "Updated atomic structure:" in line:
+                i_sym = 4
+                i_pos_min = 1 ; i_pos_max = 4
+            l += 1
+            symbols = []
+            positions = []
+            for n in range(N):
+               l += 1
+               fields = lines[l].split()
+               sym = fields[i_sym]
+               pos = [float(x) for x in fields[i_pos_min:i_pos_max]]
+               symbols.append(sym)
+               positions.append(pos)
+        elif "Total atomic forces" in line:
+            forces = []
+            for i in range(N):
+                l += 1
+                force = [float(x) for x in lines[l].split()[2:5]]
+                forces.append(force)
+        l += 1
+    
+    atoms = Atoms_with_forces(cell=cell, symbols=symbols, positions=positions)
+    atoms.forces = forces
+
+    return atoms

phonopyparser/Get_Force_Constants.py

-# Copyright 2016-2018 Fawzi Mohamed, Danio Brambila
-#
-#   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.
-
-#### phonopy parser written by Hagen-Henrik Kowalski and based on the original work of Joerg Mayer on phonopy-FHI-aims
-
-import numpy as np
-from phonopyparser.PhononModulesNomad import *
-from fnmatch import fnmatch
-import sys
-import math
-import os
-import argparse
-import logging
-
-
-from phonopy.interface.FHIaims import read_aims, write_aims, read_aims_output
-from phonopyparser.con import Control
-# Note this Phonopy is the NOMAD-lab version. Not the open source package on PYPI.
-from phonopy import Phonopy, __version__
-from phonopy.structure.symmetry import Symmetry
-from phonopy.file_IO import write_FORCE_CONSTANTS
-from phonopy.harmonic.forces import Forces
-from phonopy.harmonic.force_constants import get_force_constants
-from phonopy.units import *
-from nomadcore.unit_conversion.unit_conversion import convert_unit_function
-from nomadcore.parser_backend import *
-from nomadcore.local_meta_info import loadJsonFile, InfoKindEl
-from phonopy.structure.atoms import PhonopyAtoms as Atoms
-
-import nomad.config
-
-phonopy_version = __version__
-parser_info = {"name": "parser_phonopy", "version": "1.0"}
-
-
-class PhonopyParserWrapper():
-    """ A proper class envolop for running this parser using Noamd-FAIRD infra. """
-    def __init__(self, backend, **kwargs):
-        self.backend_factory = backend
-
-    def parse(self, mainfile):
-        logging.info('phonopy parser started')
-        logging.getLogger('nomadcore').setLevel(logging.WARNING)
-        backend = self.backend_factory("phonopy.nomadmetainfo.json")
-        # Call the old parser without a class.
-        mainDir = os.path.dirname(os.path.dirname(os.path.abspath(mainfile)))
-        cwd = os.getcwd()
-        os.chdir(mainDir)
-        try:
-            return parse_without_class(mainfile, backend)
-        finally:
-            os.chdir(cwd)
-
-
-def parse_without_class(name, backend):
-    cell_obj = read_aims("geometry.in")
-    cell = cell_obj.get_cell()
-    positions = cell_obj.get_positions()
-    symbols = np.array(cell_obj.get_chemical_symbols())
-    control = Control()
-    if (len(control.phonon["supercell"]) == 3):
-                    supercell_matrix = np.diag(control.phonon["supercell"])
-    elif (len(control.phonon["supercell"]) == 9):
-                    supercell_matrix = np.array(control.phonon["supercell"]).reshape(3,3)
-    displacement = control.phonon["displacement"]
-    sym = control.phonon["symmetry_thresh"]
-    ####
-
-    #### constructing FORCE_CONSTANTS
-    set_of_forces, phonopy_obj, Relative_Path = Collect_Forces_aims(cell_obj, supercell_matrix, displacement, sym)
-    Prep_Path = name.split("phonopy-FHI-aims-displacement-")
-    Whole_Path = []
-
-    # Try to resolve references as paths relative to the upload root.
-    try:
-        for Path in Relative_Path:
-            abs_path = "%s%s" % (Prep_Path[0], Path)
-            rel_path = abs_path.split(nomad.config.fs.staging + "/")[1].split("/", 3)[3]
-            Whole_Path.append(rel_path)
-    except Exception:
-        logging.warn("Could not resolve path to a referenced calculation within the upload.")
-    phonopy_obj.set_forces(set_of_forces)
-    phonopy_obj.produce_force_constants()
-    FC2 = phonopy_obj.get_force_constants()
-    ####
-
-    #### obtaining information about supercell
-    super_c=phonopy_obj.supercell
-    s_cell = super_c.get_cell()
-    super_pos = super_c.get_positions()
-    super_sym = np.array(super_c.get_chemical_symbols())
-    ####
-
-    #### Converting properties to Si unitis
-    converter_FC2 = convert_unit_function('eV*angstrom**-2', 'joules*meter**-2')
-    convert_angstrom = convert_unit_function('angstrom', 'meter')
-    FC2 = converter_FC2(FC2)
-    cell = convert_angstrom(cell)
-    s_cell = convert_angstrom(s_cell)
-    super_pos = convert_angstrom(super_pos)
-    positions = convert_angstrom(positions)
-    displacement = convert_angstrom(displacement)
-
-    #### parsing
-    pbc = np.array((1, 1, 1), bool)
-    Parse = backend
-    Parse.startedParsingSession(name, parser_info)
-    sRun = Parse.openSection("section_run")
-    Parse.addValue("program_name", "Phonopy")
-    Parse.addValue("program_version", phonopy_version)
-    Basesystem = Parse.openSection("section_system")
-    Parse.addArrayValues("configuration_periodic_dimensions", pbc)
-    Parse.addArrayValues("atom_labels", symbols)
-    Parse.addArrayValues("atom_positions", positions)
-    Parse.addArrayValues("simulation_cell", cell)
-    Parse.closeSection("section_system", Basesystem)
-    Supercellsystem = Parse.openSection("section_system")
-    sysrefs = Parse.openSection("section_system_to_system_refs")
-    Parse.addValue("system_to_system_kind", "subsystem")
-    Parse.addValue("system_to_system_ref", Basesystem)
-    Parse.closeSection("section_system_to_system_refs", sysrefs)
-    Parse.addArrayValues("configuration_periodic_dimensions", pbc)
-    Parse.addArrayValues("atom_labels", super_sym)
-    Parse.addArrayValues("atom_positions", super_pos)
-    Parse.addArrayValues("simulation_cell", s_cell)
-    Parse.addArrayValues("SC_matrix", supercell_matrix)
-    Parse.addValue("x_phonopy_original_system_ref", Basesystem)
-    Parse.closeSection("section_system", Supercellsystem)
-    method = Parse.openSection("section_method")
-    Parse.addValue("x_phonopy_symprec", sym)
-    Parse.addValue("x_phonopy_displacement", displacement)
-    Parse.closeSection("section_method", method)
-    results = Parse.openSection("section_single_configuration_calculation")
-    Parse.addValue("single_configuration_calculation_to_system_ref", Supercellsystem)
-    Parse.addValue("single_configuration_to_calculation_method_ref", method)
-    Parse.addArrayValues("hessian_matrix", FC2)
-    GP = Get_Properties(FC2, cell, positions, symbols, supercell_matrix, sym, displacement)
-    GP.prem_emit(Parse, results)
-    GP.prep_ref(Whole_Path, Parse)
-    Parse.closeSection("section_single_configuration_calculation", results)
-    Parse.closeSection("section_run", sRun)
-    Parse.finishedParsingSession("ParseSuccess", None)
-
-    return backend

phonopyparser/Get_Properties.py

-# Copyright 2016-2018 Fawzi Mohamed, Danio Brambila
-#
-#   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.
-
-#### phonopy parser based on the original work of Joerg Mayer on phonopy-FHI-aims
-
-import numpy as np
-from phonopyparser.PhononModulesNomad import *
-from fnmatch import fnmatch
-import sys
-import math
-import os, logging
-from phonopy.interface.FHIaims import read_aims, write_aims, read_aims_output
-from phonopyparser.con import Control
-from phonopy import Phonopy
-from phonopy.structure.symmetry import Symmetry
-from phonopy.file_IO import parse_FORCE_CONSTANTS
-from phonopy.harmonic.forces import Forces
-from phonopy.harmonic.force_constants import get_force_constants
-from phonopy.units import *
-from phonopy.structure.atoms import Atoms
-from nomadcore.unit_conversion.unit_conversion import convert_unit_function
-from nomadcore.parser_backend import *
-from nomadcore.local_meta_info import loadJsonFile, InfoKindEl
-
-parser_info = {"name": "parser_phonopy", "version": "0.1"}
-
-path = "../../../../nomad-meta-info/meta_info/nomad_meta_info/phonopy.nomadmetainfo.json"
-metaInfoPath = os.path.normpath(
-    os.path.join(os.path.dirname(os.path.abspath(__file__)), path))
-metaInfoEnv, warns = loadJsonFile(filePath=metaInfoPath,
-                                  dependencyLoader=None,
-                                  extraArgsHandling=InfoKindEl.ADD_EXTRA_ARGS,
-                                  uri=None)
-if __name__ == "__main__":
-    import sys
-    name = sys.argv[1]
-    #### Reading basic properties from JSON
-    with open(name) as FORCES:
-            data = json.load(FORCES)
-    hessian= np.array(data["sections"]["section_single_configuration_calculation-0"]["hessian_matrix"])
-    SC_matrix = np.array(data["sections"]["section_system-1"]["SC_matrix"])
-    cell = np.array(data["sections"]["section_system-0"]["simulation_cell"])
-    symbols = np.array(data["sections"]["section_system-0"]["atom_labels"])
-    positions = np.array(data["sections"]["section_system-0"]["atom_positions"])
-    displacement = np.array(data["sections"]["section_method-0"]["x_phonopy_displacement"])
-    symmetry_thresh = np.array(data["sections"]["section_method-0"]["x_phonopy_symprec"])
-    ####
-
-
-    #### omitting
-    get_properties = get_properties(hessian, cell, positions, symbols, SC_matrix, symmetry_thresh, displacement, file_name, metaInfoEnv, parser_info)
-    get_properties.omit_properties()

phonopyparser/__init__.py

@@ -12,4 +12,21 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
-from phonopyparser.Get_Force_Constants import PhonopyParserWrapper
+from phonopyparser.phonopy_calculators import PhonopyCalculatorInterface
+
+from nomad.parsing.parser import FairdiParser
+from .metainfo import m_env
+
+
+class PhonopyParser(FairdiParser):
+    def __init__(self):
+        super().__init__(
+            name='parsers/phonopy', code_name='Phonopy', code_homepage='https://phonopy.github.io/phonopy/',
+            mainfile_name_re=(r'(.*/phonopy-FHI-aims-displacement-0*1/control.in$)|(.*/phonon.yaml)')
+        )
+
+    def parse(self, filepath, archive, logger=None):
+        self._metainfo_env = m_env
+
+        interface = PhonopyCalculatorInterface(filepath, archive, logger)
+        interface.parse()

phonopyparser/__main__.py

+# 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 phonopyparser import PhonopyParser
+
+if __name__ == "__main__":
+    archive = PhonopyParser.main(sys.argv[1])
+    json.dump(archive.m_to_dict(), sys.stdout, indent=2)

phonopyparser/con.py → phonopyparser/fhiaims_io.py

@@ -12,19 +12,16 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
-#### phonopy parser based on the original work of Joerg Mayer on phonopy-FHI-aims
+# phonopy parser based on the original work of Joerg Mayer on phonopy-FHI-aims
 
 import numpy as np
-from phonopy.units import VaspToTHz as AimsToTHz, VaspToCm as AimsToCm, VaspToEv as AimsToEv, THzToCm, THzToEv
-from phonopy.interface.FHIaims import read_aims, write_aims, read_aims_output
-# from phonopy.structure.cells import get_equivalent_smallest_vectors
-from phonopy.harmonic.dynamical_matrix import  DynamicalMatrix, DynamicalMatrixNAC, get_equivalent_smallest_vectors
-from phonopy.structure.cells import get_reduced_bases
-import numpy as np
+import math
+import os
+import logging
+from fnmatch import fnmatch
 from phonopy import Phonopy
-from phonopy.structure.cells import Primitive
-
 
+from phonopyparser.FHIaims import read_aims_output
 
 class Control:
     def __init__(self, file=None):
@@ -82,17 +79,89 @@ class Control:
                             delta = (0.1, 0.1, 0.1)
                         if (delta[0] == 0.0) and (delta[1] == 0.0) and (delta[2] == 0.0):
                             raise Exception("evaluation of frequencies with non-analytic corrections must be shifted by delta away from Gamma")
-                        parameters = { "file" : fields[2],
-                                       "method" : fields[3].lower(),
-                                       "delta" : delta }
+                        parameters = {
+                            "file": fields[2], "method": fields[3].lower(), "delta": delta}
                         self.phonon["nac"].update(parameters)
 
-        except Exception:
-            #print (line,)
-            #print ("|-> line triggered exception: ") + str(Exception)
-            raise
+        except Exception as e:
+            raise(e)
+
         # supercell is mandatory for all what follows
         if not self.phonon["supercell"]:
             raise Exception("no supercell specified in %s" % self.file)
         f.close()
 
+
+def clean_position(scaled_positions):
+    scaled_positions = list(scaled_positions)
+    for sp in range(len(scaled_positions)):
+        for i in range(len(scaled_positions[sp])):
+            if np.float(np.round(scaled_positions[sp][i], 7)) >= 1:
+                scaled_positions[sp][i] -= 1.0
+            elif scaled_positions[sp][i] <= -1e-5:
+                scaled_positions[sp][i] += 1.0
+    scaled_positions = np.array(scaled_positions)
+
+    return scaled_positions
+
+
+def read_forces_aims(cell_obj, supercell_matrix, displacement, sym, tol=1e-6, logger=None):
+    if logger is None:
+        logger = logging
+
+    phonopy_obj = Phonopy(cell_obj, supercell_matrix, symprec=sym)
+    phonopy_obj.generate_displacements(distance=displacement)
+    supercells = phonopy_obj.get_supercells_with_displacements()
+    directories = []
+    digits = int(math.ceil(math.log(len(supercells) + 1, 10))) + 1
+    for i in range(len(supercells)):
+        directories.append(("phonopy-FHI-aims-displacement-%0" + str(digits) + "d") % (i + 1))
+    set_of_forces = []
+    Relative_Path = []
+    for directory, supercell in zip(directories, supercells):
+        aims_out = os.path.join(directory, directory + ".out")
+        if not os.path.isfile(aims_out):
+            logger.warn("!!! file not found: %s" % aims_out)
+            cwd = os.getcwd()
+            con_list = os.listdir(cwd)
+            check_var = False
+            for name in con_list:
+                if fnmatch(name, '*.out'):
+                    aims_out = '%s/%s' % (directory, name)
+                    logger.warn(
+                        "Your file seems to have a wrong name proceeding with %s" % aims_out
+                    )
+                    check_var = True
+                    break
+            if not check_var:
+                raise Exception("No phonon calculations found")
+            os.chdir("../")
+        Relative_Path.append(aims_out)
+        supercell_calculated = read_aims_output(aims_out)
+        if (
+            (supercell_calculated.get_number_of_atoms() == supercell.get_number_of_atoms()) and
+            (supercell_calculated.get_atomic_numbers() == supercell.get_atomic_numbers()).all() and
+            (abs(supercell_calculated.get_positions() - supercell.get_positions()) < tol).all() and
+            (abs(supercell_calculated.get_cell() - supercell.get_cell()) < tol).all()):
+            # read_aims_output reads in forces from FHI-aims output as list structure,
+            # but further processing below requires numpy array
+            forces = np.array(supercell_calculated.get_forces())
+            drift_force = forces.sum(axis=0)
+            for force in forces:
+                force -= drift_force / forces.shape[0]
+            set_of_forces.append(forces)
+        elif (
+            (supercell_calculated.get_number_of_atoms() == supercell.get_number_of_atoms()) and
+            (supercell_calculated.get_atomic_numbers() == supercell.get_atomic_numbers()).all() and
+            (abs(clean_position(supercell_calculated.get_scaled_positions()) - clean_position(supercell.get_scaled_positions())) < tol).all() and
+            (abs(supercell_calculated.get_cell() - supercell.get_cell()) < tol).all()):
+            logger.warn("!!! there seems to be a rounding error")
+            forces = np.array(supercell_calculated.get_forces())
+            drift_force = forces.sum(axis=0)
+            for force in forces:
+                force -= drift_force / forces.shape[0]
+            set_of_forces.append(forces)
+        else:
+            raise Exception("calculated varies from expected supercell in FHI-aims output %s" % aims_out)
+
+    return set_of_forces, phonopy_obj, Relative_Path

phonopyparser/phonopy_calculators.py

+import os
+import numpy as np
+import phonopy
+from phonopy.units import THzToEv
+
+from phonopyparser.fhiaims_io import Control, read_forces_aims
+from phonopyparser.phonopy_properties import PhononProperties
+from phonopyparser.FHIaims import read_aims
+
+from nomadcore.unit_conversion.unit_conversion import convert_unit_function
+
+from nomad.datamodel.metainfo.public import section_run, section_system,\
+    section_system_to_system_refs, section_method, section_single_configuration_calculation,\
+    section_k_band, section_k_band_segment, section_dos, section_frame_sequence,\
+    section_thermodynamical_properties, section_sampling_method, Workflow, Phonon, \
+    section_calculation_to_calculation_refs
+
+
+class PhonopyCalculatorInterface:
+    def __init__(self, filepath, archive, logger=None):
+        self.filepath = os.path.abspath(filepath)
+        self.archive = archive
+        self.logger = logger
+        self._phonopy_obj = None
+        self.calculator = None
+        if 'control.in' in filepath:
+            self.calculator = 'fhi-aims'
+        elif 'phonon.yaml' in filepath:
+            self.calculator = 'vasp'
+        self.references = []
+
+    @property
+    def phonopy_obj(self):
+        if self._phonopy_obj is None:
+            if self.calculator == 'fhi-aims':
+                self._build_phonopy_object_fhi_aims()
+            elif self.calculator == 'vasp':
+                self._build_phonopy_object_vasp()
+        return self._phonopy_obj
+
+    def _build_phonopy_object_vasp(self):
+        cwd = os.getcwd()
+        os.chdir(os.path.dirname(self.filepath))
+
+        phonopy_obj = phonopy.load('phonon.yaml')
+        os.chdir(cwd)
+
+        self._phonopy_obj = phonopy_obj
+
+    def _build_phonopy_object_fhi_aims(self):
+        cwd = os.getcwd()
+        os.chdir(os.path.dirname(os.path.dirname(self.filepath)))
+
+        cell_obj = read_aims("geometry.in")
+        control = Control()
+        if (len(control.phonon["supercell"]) == 3):
+            supercell_matrix = np.diag(control.phonon["supercell"])
+        elif (len(control.phonon["supercell"]) == 9):
+            supercell_matrix = np.array(control.phonon["supercell"]).reshape(3, 3)
+        displacement = control.phonon["displacement"]
+        sym = control.phonon["symmetry_thresh"]
+
+        set_of_forces, phonopy_obj, Relative_Path = read_forces_aims(
+            cell_obj, supercell_matrix, displacement, sym)
+        # Prep_Path = name.split("phonopy-FHI-aims-displacement-")
+
+        phonopy_obj.set_forces(set_of_forces)
+        phonopy_obj.produce_force_constants()
+
+        os.chdir(cwd)
+
+        self._phonopy_obj = phonopy_obj
+
+    def parse_bandstructure(self):
+        freqs, bands, bands_labels = self.properties.get_bandstructure()
+
+        # convert THz to eV
+        freqs = freqs * THzToEv
+
+        # convert eV to J
+        eVtoJoules = convert_unit_function('eV', 'joules')
+        freqs = eVtoJoules(freqs)
+
+        sec_scc = self.archive.section_run[0].section_single_configuration_calculation[0]
+
+        sec_k_band = sec_scc.m_create(section_k_band)
+        sec_k_band.band_structure_kind = 'vibrational'
+
+        for i in range(len(freqs)):
+            freq = np.expand_dims(freqs[i], axis=0)
+            sec_k_band_segment = sec_k_band.m_create(section_k_band_segment)
+            sec_k_band_segment.band_energies = freq
+            sec_k_band_segment.band_k_points = bands[i]
+            sec_k_band_segment.band_segm_labels = bands_labels[i]
+
+    def parse_dos(self):
+        f, dos = self.properties.get_dos()
+
+        # To match the shape given in meta data another dimension is added to the
+        # array (spin degress of fredom is 1)
+        dos = np.expand_dims(dos, axis=0)
+
+        # convert THz to eV to Joules
+        eVtoJoules = convert_unit_function('eV', 'joules')
+        f = f * THzToEv
+        f = eVtoJoules(f)
+
+        sec_scc = self.archive.section_run[0].section_single_configuration_calculation[0]
+        sec_dos = sec_scc.m_create(section_dos)
+        sec_dos.dos_kind = 'vibrational'
+        sec_dos.dos_values = dos
+        sec_dos.dos_energies = f
+
+    def parser_thermodynamical_properties(self):
+        T, fe, _, cv = self.properties.get_thermodynamical_properties()
+
+        n_atoms = len(self.phonopy_obj.unitcell)
+        n_atoms_supercell = len(self.phonopy_obj.supercell)
+
+        fe = fe / n_atoms
+
+        # The thermodynamic properties are reported by phonopy for the base
+        # system. Since the values in the metainfo are stored per the referenced
+        # system, we need to multiple by the size factor between the base system
+        # and the supersystem used in the calculations.
+        cv = cv * (n_atoms_supercell / n_atoms)
+
+        # convert to SI units
+        eVtoJoules = convert_unit_function('eV', 'joules')
+        eVperKtoJoules = convert_unit_function('eV*K**-1', 'joules*K**-1')
+        fe = eVtoJoules(fe)
+        cv = eVperKtoJoules(cv)
+
+        sec_run = self.archive.section_run[0]
+        sec_scc = sec_run.section_single_configuration_calculation
+
+        sec_frame_sequence = sec_run.m_create(section_frame_sequence)
+        sec_frame_sequence.frame_sequence_local_frames_ref = sec_scc
+
+        sec_thermo_prop = sec_frame_sequence.m_create(section_thermodynamical_properties)
+        sec_thermo_prop.thermodynamical_property_temperature = T
+        sec_thermo_prop.vibrational_free_energy_at_constant_volume = fe
+        sec_thermo_prop.thermodynamical_property_heat_capacity_C_v = cv
+
+        sec_sampling_method = sec_run.m_create(section_sampling_method)
+        sec_sampling_method.sampling_method = 'taylor_expansion'
+        sec_sampling_method.sampling_method_expansion_order = 2
+
+        sec_frame_sequence.frame_sequence_to_sampling_ref = sec_sampling_method
+
+    def parse_ref(self):
+        sec_scc = self.archive.section_run[0].section_single_configuration_calculation[0]
+        sec_calc_refs = sec_scc.m_create(section_calculation_to_calculation_refs)
+        sec_calc_refs.calculation_to_calculation_kind = 'source_calculation'
+        for ref in self.references:
+            sec_calc_refs.calculation_to_calculation_external_url = ref
+
+    def parse(self):
+        phonopy_obj = self.phonopy_obj
+        self.properties = PhononProperties(self.phonopy_obj, self.logger)
+
+        pbc = np.array((1, 1, 1), bool)
+
+        unit_cell = self.phonopy_obj.unitcell.get_cell()
+        unit_pos = self.phonopy_obj.unitcell.get_positions()
+        unit_sym = np.array(self.phonopy_obj.unitcell.get_chemical_symbols())
+
+        super_cell = self.phonopy_obj.supercell.get_cell()
+        super_pos = self.phonopy_obj.supercell.get_positions()
+        super_sym = np.array(self.phonopy_obj.supercell.get_chemical_symbols())
+
+        # convert to SI
+        convert_fc = convert_unit_function('eV*angstrom**-2', 'joules*meter**-2')
+        convert_angstrom = convert_unit_function('angstrom', 'meter')
+
+        unit_cell = convert_angstrom(unit_cell)
+        unit_pos = convert_angstrom(unit_pos)
+
+        super_cell = convert_angstrom(super_cell)
+        super_pos = convert_angstrom(super_pos)
+
+        displacement = np.linalg.norm(phonopy_obj.displacements[0][1:])
+        displacement = convert_angstrom(displacement)
+        supercell_matrix = phonopy_obj.supercell_matrix
+        sym_tol = phonopy_obj.symmetry.tolerance
+        force_constants = phonopy_obj.get_force_constants()
+        force_constants = convert_fc(force_constants)
+
+        sec_run = self.archive.m_create(section_run)
+        sec_run.program_name = 'Phonopy'
+        sec_run.program_version = phonopy.__version__
+        sec_system_unit = sec_run.m_create(section_system)
+        sec_system_unit.configuration_periodic_dimensions = pbc
+        sec_system_unit.atom_labels = unit_sym
+        sec_system_unit.atom_positions = unit_pos
+        sec_system_unit.simulation_cell = unit_cell
+
+        sec_system = sec_run.m_create(section_system)
+        sec_system_to_system_refs = sec_system.m_create(section_system_to_system_refs)
+        sec_system_to_system_refs.system_to_system_kind = 'subsystem'
+        sec_system_to_system_refs.system_to_system_ref = sec_system_unit
+        sec_system.configuration_periodic_dimensions = pbc
+        sec_system.atom_labels = super_sym
+        sec_system.atom_positions = super_pos
+        sec_system.simulation_cell = super_cell
+        sec_system.SC_matrix = supercell_matrix
+        sec_system.x_phonopy_original_system_ref = sec_system_unit
+
+        sec_method = sec_run.m_create(section_method)
+        sec_method.x_phonopy_symprec = sym_tol
+        sec_method.x_phonopy_displacement = displacement
+
+        sec_scc = sec_run.m_create(section_single_configuration_calculation)
+        sec_scc.single_configuration_calculation_to_system_ref = sec_system
+        sec_scc.single_configuration_to_calculation_method_ref = sec_method
+        sec_scc.hessian_matrix = force_constants
+
+        self.parse_bandstructure()
+        self.parse_dos()
+        self.parser_thermodynamical_properties()
+        self.parse_ref()
+
+        sec_workflow = self.archive.m_create(Workflow)
+        sec_workflow.workflow_type = 'phonon'
+        sec_phonon = sec_workflow.m_create(Phonon)
+        sec_phonon.force_calculator = self.calculator
+        vol = np.dot(unit_cell[0], np.cross(unit_cell[1], unit_cell[2]))
+        sec_phonon.mesh_density = np.prod(self.properties.mesh) / vol
+        n_imaginary = np.count_nonzero(self.properties.frequencies < 0)
+        sec_phonon.n_imaginary_frequencies = n_imaginary
+        if phonopy_obj.nac_params:
+            sec_phonon.with_non_analytic_correction = True

phonopyparser/phonopy_properties.py

+import numpy as np
+import logging
+
+from ase.geometry import crystal_structure_from_cell
+from ase.dft.kpoints import special_paths, parse_path_string
+try:
+    from ase.dft.kpoints import special_points
+except ImportError:
+    from ase.dft.kpoints import sc_special_points as special_points
+
+from phonopy.phonon.band_structure import BandStructure
+from phonopy.units import EvTokJmol, VaspToTHz
+
+
+def generate_kpath_ase(cell, symprec):
+
+    eig_val_max = np.real(np.linalg.eigvals(cell)).max()
+    eps = eig_val_max * symprec
+
+    lattice = crystal_structure_from_cell(cell, eps)
+    paths = special_paths.get(lattice, None)
+    if paths is None:
+        paths = special_paths['orthorhombic']
+    paths = parse_path_string(special_paths[lattice])
+    points = special_points[lattice]
+    k_points = []
+    for p in paths:
+        k_points.append([points[k] for k in p])
+        for index in range(len(p)):
+            if p[index] == 'G':
+                p[index] = 'Γ'
+    parameters = []
+    for h, seg in enumerate(k_points):
+        for i, path in enumerate(seg):
+            parameter = {}
+            parameter['npoints'] = 100
+            parameter['startname'] = paths[h][i]
+            if i == 0 and len(seg) > 2:
+                parameter['kstart'] = path
+                parameter['kend'] = seg[i + 1]
+                parameter['endname'] = paths[h][i + 1]
+                parameters.append(parameter)
+            elif i == (len(seg) - 2):
+                parameter['kstart'] = path
+                parameter['kend'] = seg[i + 1]
+                parameter['endname'] = paths[h][i + 1]
+                parameters.append(parameter)
+                break
+            else:
+                parameter['kstart'] = path
+                parameter['kend'] = seg[i + 1]
+                parameter['endname'] = paths[h][i + 1]
+                parameters.append(parameter)
+    return parameters
+
+
+class PhononProperties():
+    def __init__(self, phonopy_obj, logger=None):
+        if logger is None:
+            logger = logging
+
+        self.phonopy_obj = phonopy_obj
+
+        self.n_atoms = len(phonopy_obj.unitcell)
+
+        mesh_density = (2 * 80 ** 3) / self.n_atoms
+        mesh_number = np.round(mesh_density**(1. / 3.))
+        self.mesh = [mesh_number, mesh_number, mesh_number]
+
+        self.n_atoms_supercell = len(phonopy_obj.supercell)
+
+    def get_bandstructure(self):
+        phonopy_obj = self.phonopy_obj
+
+        frequency_unit_factor = VaspToTHz
+        is_eigenvectors = False
+
+        unit_cell = phonopy_obj.unitcell.get_cell()
+        sym_tol = phonopy_obj.symmetry.tolerance
+        parameters = generate_kpath_ase(unit_cell, sym_tol)
+
+        # Distances calculated in phonopy.band_structure.BandStructure object
+        # are based on absolute positions of q-points in reciprocal space
+        # as calculated by using the cell which is handed over during instantiation.
+        # Fooling that object by handing over a "unit cell" diag(1,1,1) instead clashes
+        # with calculation of non-analytical terms.
+        # Hence generate appropriate distances and special k-points list based on fractional
+        # coordinates in reciprocal space (to keep backwards compatibility with previous
+        # FHI-aims phonon implementation).
+        bands = []
+        bands_distances = []
+        distance = 0.0
+        bands_special_points = [distance]
+        bands_labels = []
+        label = parameters[0]["startname"]
+        for b in parameters:
+            kstart = np.array(b["kstart"])
+            kend = np.array(b["kend"])
+            npoints = b["npoints"]
+            dk = (kend - kstart) / (npoints - 1)
+            bands.append([(kstart + dk * n) for n in range(npoints)])
+            dk_length = np.linalg.norm(dk)
+
+            for n in range(npoints):
+                bands_distances.append(distance + dk_length * n)
+
+            distance += dk_length * (npoints - 1)
+            bands_special_points.append(distance)
+            label = [b["startname"], b["endname"]]
+            bands_labels.append(label)
+
+        bs_obj = BandStructure(
+            bands, phonopy_obj.dynamical_matrix, with_eigenvectors=is_eigenvectors,
+            factor=frequency_unit_factor)
+
+        freqs = bs_obj.get_frequencies()
+
+        return np.array(freqs), np.array(bands), np.array(bands_labels)
+
+    def get_dos(self):
+        phonopy_obj = self.phonopy_obj
+        mesh = self.mesh
+
+        phonopy_obj.set_mesh(mesh, is_gamma_center=True)
+        q_points = phonopy_obj.get_mesh()[0]
+        phonopy_obj.set_qpoints_phonon(q_points, is_eigenvectors=False)
+
+        frequencies = phonopy_obj.get_qpoints_phonon()[0]
+        self.frequencies = np.array(frequencies)
+        min_freq = min(np.ravel(frequencies))
+        max_freq = max(np.ravel(frequencies)) + max(np.ravel(frequencies)) * 0.05
+
+        phonopy_obj.set_total_DOS(
+            freq_min=min_freq, freq_max=max_freq, tetrahedron_method=True)
+        f, dos = phonopy_obj.get_total_DOS()
+
+        return f, dos
+
+    def get_thermodynamical_properties(self):
+        phonopy_obj = self.phonopy_obj
+
+        t_max = 100
+        t_min = 0
+        t_step = 10
+        mesh = self.mesh
+
+        phonopy_obj.set_mesh(mesh, is_gamma_center=True)
+        phonopy_obj.set_thermal_properties(t_step=t_step, t_max=t_max, t_min=t_min)
+        T, fe, entropy, cv = phonopy_obj.get_thermal_properties()
+        kJmolToEv = 1.0 / EvTokJmol
+        fe = fe * kJmolToEv
+        JmolToEv = kJmolToEv / 1000
+        cv = JmolToEv * cv
+        return T, fe, entropy, cv

requirements.txt

-phonopy@git+https://gitlab.mpcdf.mpg.de/nomad-lab/phonopy@master#egg=phonopy
\ No newline at end of file
+phonopy

setup.py

@@ -25,7 +25,7 @@ def main():
         packages=find_packages(),
         install_requires=[
             'nomadcore',
-            'phonopy-nomad-lab'
+            'phonopy'
         ],
     )