Enabled the reading of system geometries from input file if they are not...

Enabled the reading of system geometries from input file if they are not present in the output file.

parser/parser-cp2k/cp2kparser/generic/inputparsing.py

@@ -53,7 +53,7 @@ class CP2KInput(object):
         else:
             parameter.value = value
 
-    def set_keyword(self, path, value):
+    def set_keyword(self, path, value, full):
         keyword, section = self.get_keyword_and_section(path)
         # If keyword found, put data in there
         if keyword and section:
@@ -63,8 +63,7 @@ class CP2KInput(object):
             split_path = path.rsplit("/", 1)
             keyword = split_path[1]
             if section.default_keyword is not None:
-                # print "Saving default keyword at path '{}'".format(path)
-                section.default_keyword.value += keyword + " " + value + "\n"
+                section.default_keyword.value += full + "\n"
             else:
                 message = "The CP2K input does not contain the keyword {}, and there is no default keyword for the section {}".format(path, split_path[0])
                 logger.warning(message)
@@ -140,6 +139,7 @@ class CP2KInput(object):
             return (None, None)
         if section.section_parameter is not None:
             parameter = section.section_parameter
+            parameter.accessed = section.accessed
             return (parameter, section)
         else:
             return (None, section)
@@ -337,12 +337,15 @@ class Keyword(InputObject):
 
     def get_unit(self):
 
-        # Decode the unit and the value if not done before
-        if self.default_unit:
-            if not self.unit:
+        if self.unit is None:
             self.decode_cp2k_unit_and_value()
+            if self.unit is not None:
                 return self.unit
+            elif self.default_unit is not None:
+                return self.default_unit
         else:
+            return self.unit
+
         logger.error("The keyword '{}' does not have a unit.".format(self.default_name))
         return None
 
@@ -371,12 +374,26 @@ class SectionParameters(InputObject):
     Section parameters are the short values that can be added right after a
     section name, e.g. &PRINT ON, where ON is the section parameter.
     """
-    __slots__ = ['lone_keyword_value']
+    __slots__ = ['lone_keyword_value', 'accessed']
 
     def __init__(self, default_value, lone_keyword_value):
         super(SectionParameters, self).__init__("SECTION_PARAMETERS")
         self.default_value = default_value
         self.lone_keyword_value = lone_keyword_value
+        self.accessed = None
+
+    def get_value(self):
+        """Returns the value for this section parameter. Uses the user given
+        value primarily, if the section is not used, return the default value
+        and if the section is defined but without explicit section parameter
+        returns the lone keyword value.
+        """
+        if self.value is not None:
+            return self.value
+        elif self.accessed:
+            return self.lone_keyword_value
+        else:
+            return self.default_value
 
 
 class DefaultKeyword(InputObject):

parser/parser-cp2k/cp2kparser/versions/cp2k262/inputparser.py

@@ -7,6 +7,7 @@ import logging
 import pickle
 import numpy as np
 from nomadcore.baseclasses import AbstractBaseParser
+import nomadcore.configurationreading
 from cp2kparser.generic.inputparsing import metainfo_data_prefix, metainfo_section_prefix
 from pint import UnitRegistry
 import ase
@@ -89,14 +90,14 @@ class CP2KInputParser(AbstractBaseParser):
         print_level = self.input_tree.get_keyword("GLOBAL/PRINT_LEVEL")
         multiple_unit_cell = self.input_tree.get_keyword("FORCE_EVAL/SUBSYS/MULTIPLE_UNIT_CELL")
         top_multiple_unit_cell = self.input_tree.get_keyword("FORCE_EVAL/SUBSYS/TOPOLOGY/MULTIPLE_UNIT_CELL")
-        centering = self.input_tree.get_section("FORCE_EVAL/SUBSYS/TOPOLOGY/CENTER_COORDINATES")
+        centering = self.input_tree.get_section("FORCE_EVAL/SUBSYS/TOPOLOGY/CENTER_COORDINATES").section_parameter.get_value()
+        center_point = self.input_tree.get_keyword("FORCE_EVAL/SUBSYS/TOPOLOGY/CENTER_COORDINATES/CENTER_POINT")
         if multiple_unit_cell is None:
             multiple_unit_cell = np.array([1, 1, 1])
 
         if print_level == "LOW" \
            and np.array_equal(multiple_unit_cell, [1, 1, 1]) \
-           and np.array_equal(top_multiple_unit_cell, [1, 1, 1]) \
-           and not centering.accessed:
+           and np.array_equal(top_multiple_unit_cell, [1, 1, 1]):
 
             cell_section = self.input_tree.get_section("FORCE_EVAL/SUBSYS/CELL")
             a_obj = cell_section.get_keyword_object("A")
@@ -133,11 +134,16 @@ class CP2KInputParser(AbstractBaseParser):
                 self.cache_service["simulation_cell"] = cell
                 self.cache_service["lattice_vectors"] = cell
 
+            # Parse the coordinates from input file or xyz file
             coord = self.input_tree.get_section("FORCE_EVAL/SUBSYS/COORD")
+            coord_filename = self.input_tree.get_keyword("FORCE_EVAL/SUBSYS/TOPOLOGY/COORD_FILE_NAME")
+            coord_format = self.input_tree.get_keyword("FORCE_EVAL/SUBSYS/TOPOLOGY/COORD_FILE_FORMAT")
             scaled = coord.get_keyword("SCALED")
             coords = coord.default_keyword.value
-            print(coords)
             coord_unit = coord.get_keyword("UNIT")
+
+            final_pos = None
+            final_lab = None
             if len(coords) != 0:
                 lines = coords.split("\n")
                 pos = []
@@ -152,20 +158,53 @@ class CP2KInputParser(AbstractBaseParser):
                         labels.append(i_lab)
                         pos.append([i_x, i_y, i_z])
 
-                pos = np.array(pos)
+                final_pos = np.array(pos)
+                final_lab = np.array(labels)
+
+            elif coord_filename is not None:
+
+                extension_map = {
+                    "XYZ": "xyz",
+                }
+                coord_extension = extension_map.get(coord_format)
+                if coord_extension is not None:
+                    abs_filename = self.file_service.get_absolute_path_to_file(coord_filename)
+                    try:
+                        atoms = ase.io.read(abs_filename, format=coord_extension)
+                        final_pos = atoms.get_positions()
+                        final_lab = atoms.get_chemical_symbols()
+                    except StopIteration:
+                        logger.error("Error in reading the structure file specified in the input file.")
+                else:
+                    logger.error(
+                        "Could not read the structure file specified in the "
+                        "input file because the file format is not supported."
+                    )
+
+            # Make the positions cartesian if they are scaled
+            if final_pos is not None and final_lab is not None:
+                if scaled == "T" or scaled == "TRUE":
+                    atoms = ase.Atoms(
+                        scaled_positions=final_pos,
+                        symbols=final_lab,
+                        cell=cell,
+                    )
+                    final_pos = atoms.get_positions()
+
+                # Center the atoms if requested.
+                if centering == "T" or centering == "TRUE":
+                    atoms = ase.Atoms(
+                        positions=final_pos,
+                        symbols=final_lab,
+                        cell=cell,
+                    )
+                    atoms.center(about=center_point)
+                    final_pos = atoms.get_positions()
+
                 coord_unit = self.get_pint_unit_string(coord_unit)
-                pos = (pos * ureg(coord_unit)).to("angstrom").magnitude
-                labels = np.array(labels)
-                if scaled == "F" or scaled == "FALSE":
-                    self.cache_service["atom_positions"] = pos
-                    self.cache_service["atom_labels"] = labels
-
-            # self.cache_service["simulation_cell"] = cell
-            # self.cache_service["lattice_vectors"] = cell
-            # coordinates = self.input_tree.get_section("FORCE_EVAL/SUBSYS/COORD")
-            # print(coordinates.get_section_parameter())
-            # self.cache_service["atom_positions"] = coordinates
-            # self.cache_service["atom_labels"] = labels
+                final_pos = (final_pos * ureg(coord_unit)).to("angstrom").magnitude
+                self.cache_service["atom_positions"] = final_pos
+                self.cache_service["atom_labels"] = final_lab
 
         # Parse the used XC_functionals and their parameters
         xc = self.input_tree.get_section("FORCE_EVAL/DFT/XC/XC_FUNCTIONAL")
@@ -302,8 +341,8 @@ class CP2KInputParser(AbstractBaseParser):
         self.cache_service["velocity_unit"] = pint_vel_unit
 
         #=======================================================================
+        # OLD: This information is retrieved from output file
         # See if some more exotic calculation is requested (e.g. MP2, DFT+U, GW, RPA)
-
         # Search for a WF_CORRELATION section
         # correlation = self.input_tree.get_section("FORCE_EVAL/DFT/XC/WF_CORRELATION")
         # method = "DFT"
@@ -470,14 +509,12 @@ class CP2KInputParser(AbstractBaseParser):
             # Contents (keywords, default keywords)
             else:
                 split = line.split(None, 1)
-                print(split)
                 if len(split) <= 1:
                     keyword_value = ""
                 else:
                     keyword_value = split[1]
                 keyword_name = split[0].upper()
-                # print(keyword_value)
-                self.input_tree.set_keyword(path + "/" + keyword_name, keyword_value)
+                self.input_tree.set_keyword(path + "/" + keyword_name, keyword_value, line)
 
     def fill_metadata(self):
         """Goes through the input data and pushes everything to the

regtests/regtests.py

@@ -251,26 +251,47 @@ class TestConfigurationPeriodicDimensions(unittest.TestCase):
 class TestGeometryInfo(unittest.TestCase):
     """Tests for a CP2K calculation with RUN_TYPE ENERGY_FORCE.
     """
-    # def test_input_file_vectors(self):
-        # """Test geometry given in input file only
-        # """
-        # results = get_result("geometries/input_vectors")
-        # simulation_cell = results["simulation_cell"]
-        # lattice_vectors = results["lattice_vectors"]
-
-        # expected = np.array([
-            # [5.430697500, 0.000000000, 0.000000000],
-            # [0.000000000, 5.430697500, 0.000000000],
-            # [0.000000000, 0.000000000, 5.430697500]
-        # ])*1e-10
-
-        # self.assertTrue(np.allclose(simulation_cell, expected, atol=1e-15, rtol=0))
-        # self.assertTrue(np.allclose(lattice_vectors, expected, atol=1e-15, rtol=0))
+    def test_input_file_vectors(self):
+        """Test geometry given in input file only
+        """
+        results = get_result("geometries/input_vectors")
+
+        # Test the cell information
+        simulation_cell = results["simulation_cell"]
+        lattice_vectors = results["lattice_vectors"]
+
+        expected = np.array([
+            [5.430697500, 0.000000000, 0.000000000],
+            [0.000000000, 5.430697500, 0.000000000],
+            [0.000000000, 0.000000000, 5.430697500]
+        ])*1e-10
+
+        self.assertTrue(np.allclose(simulation_cell, expected, atol=1e-15, rtol=0))
+        self.assertTrue(np.allclose(lattice_vectors, expected, atol=1e-15, rtol=0))
+
+        # Test the coordinate and label information
+        pos = results["atom_positions"]
+        labels = results["atom_labels"]
+        expected_pos = np.array([
+            [0.000000000, 0.000000000, 0.000000000],
+            [0.000000000, 2.715348700, 2.715348700],
+            [2.715348700, 2.715348700, 0.000000000],
+            [2.715348700, 0.000000000, 2.715348700],
+            [4.073023100, 1.357674400, 4.073023100],
+            [1.357674400, 1.357674400, 1.357674400],
+            [1.357674400, 4.073023100, 4.073023100],
+            [4.073023100, 4.073023100, 1.357674400]
+        ])*1e-10
+        self.assertTrue(np.allclose(pos, expected_pos, atol=1e-15, rtol=0))
+        expected_lab = np.array(["Si", "Si", "Si", "Si", "Si", "Si", "Si", "Si"])
+        self.assertTrue(np.array_equal(labels, expected_lab))
 
     def test_input_file_parameters(self):
         """Test geometry given as angles and lengths and with radians as unit.
         """
         results = get_result("geometries/input_parameters")
+
+        # Test the cell information
         simulation_cell = results["simulation_cell"]
         lattice_vectors = results["lattice_vectors"]
 
@@ -283,10 +304,58 @@ class TestGeometryInfo(unittest.TestCase):
         self.assertTrue(np.allclose(simulation_cell, expected, atol=1e-15, rtol=0))
         self.assertTrue(np.allclose(lattice_vectors, expected, atol=1e-15, rtol=0))
 
+        # Test the coordinate and label information
         pos = results["atom_positions"]
         labels = results["atom_labels"]
-        # print(pos)
-        # print(labels)
+        expected_pos = np.array([
+            [0.000000000, 0.000000000, 0.000000000],
+            [0.000000000, 2.715348700, 2.715348700],
+            [2.715348700, 2.715348700, 0.000000000],
+            [2.715348700, 0.000000000, 2.715348700],
+            [4.073023100, 1.357674400, 4.073023100],
+            [1.357674400, 1.357674400, 1.357674400],
+            [1.357674400, 4.073023100, 4.073023100],
+            [4.073023100, 4.073023100, 1.357674400]
+        ])*1e-10
+        self.assertTrue(np.allclose(pos, expected_pos, atol=1e-15, rtol=0))
+        expected_lab = np.array(["Si", "Si", "Si", "Si", "Si", "Si", "Si", "Si"])
+        self.assertTrue(np.array_equal(labels, expected_lab))
+
+    def test_input_file_xyz_centered(self):
+        """Test geometry file given as xyz and centered.
+        """
+        results = get_result("geometries/input_xyz_centered")
+
+        # Test the cell information
+        simulation_cell = results["simulation_cell"]
+        lattice_vectors = results["lattice_vectors"]
+
+        expected = np.array([
+            [8, 0, 0],
+            [0, 8, 0],
+            [0, 0, 8]
+        ])*1e-10
+        self.assertTrue(np.allclose(simulation_cell, expected, atol=1e-15, rtol=0))
+        self.assertTrue(np.allclose(lattice_vectors, expected, atol=1e-15, rtol=0))
+
+        # Test the coordinate and label information
+        pos = results["atom_positions"]
+        labels = results["atom_labels"]
+        expected_pos = np.array([
+            [4.000000, 4.000000, 4.500000],
+            [4.000000, 4.000000, 3.500000],
+        ])*1e-10
+        self.assertTrue(np.allclose(pos, expected_pos, atol=1e-15, rtol=0))
+        expected_lab = np.array(["H", "H"])
+        self.assertTrue(np.array_equal(labels, expected_lab))
+
+        # Test that other info is OK
+        program_name = results["program_name"]
+        esm = results["electronic_structure_method"]
+        basis_set = results["program_basis_set_type"]
+        self.assertEqual(program_name, "CP2K")
+        self.assertEqual(esm, "DFT")
+        self.assertEqual(basis_set, "gaussians")
 
 
 class TestEnergyForce(unittest.TestCase):
@@ -1192,27 +1261,27 @@ if __name__ == '__main__':
 
     for VERSION in VERSIONS:
         suites = []
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestErrors))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestUnknownInput))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestXCFunctional))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestEnergyForce))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestErrors))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestUnknownInput))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestXCFunctional))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestEnergyForce))
         suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeometryInfo))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestStressTensorMethods))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSelfInteractionCorrectionMethod))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestConfigurationPeriodicDimensions))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSCFConvergence))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestForceFiles))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestPreprocessor))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOpt))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOptTrajFormats))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOptOptimizers))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOptTrajectory))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestRestart))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMD))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMDPrintSettings))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMDEnsembles))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMDEnsembles))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestVDWMethod))
-        # suites.append(unittest.TestLoader().loadTestsFromTestCase(TestElectronicStructureMethod))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestStressTensorMethods))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSelfInteractionCorrectionMethod))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestConfigurationPeriodicDimensions))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSCFConvergence))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestForceFiles))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestPreprocessor))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOpt))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOptTrajFormats))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOptOptimizers))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestGeoOptTrajectory))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestRestart))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMD))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMDPrintSettings))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMDEnsembles))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestMDEnsembles))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestVDWMethod))
+        suites.append(unittest.TestLoader().loadTestsFromTestCase(TestElectronicStructureMethod))
         alltests = unittest.TestSuite(suites)
         unittest.TextTestRunner(verbosity=0).run(alltests)