From fe0d9f5c34f1d40f5a99257de582283a37022a44 Mon Sep 17 00:00:00 2001
From: "Himanen, Lauri (himanel1)" <lauri.himanen@aalto.fi>
Date: Wed, 18 May 2016 21:39:50 +0300
Subject: [PATCH] Added more testing tools, updated stress tensor parsing.

---
 .gitignore                                    |  1 +
 .../versions/cp2k262/commonmatcher.py         | 87 +++++++++++--------
 .../versions/cp2k262/inputparser.py           | 14 +++
 .../versions/cp2k262/singlepointparser.py     | 42 +++++++--
 test/unittests/cp2k_2.6.2/run_tests.py        | 60 ++++++++++++-
 5 files changed, 157 insertions(+), 47 deletions(-)

diff --git a/.gitignore b/.gitignore
index 85d5d27..d36cb68 100644
--- a/.gitignore
+++ b/.gitignore
@@ -64,6 +64,7 @@ test/**/*.restart.bak-3
 test/**/*.wfn.bak-1
 test/**/*.wfn.bak-2
 test/**/*.wfn.bak-3
+test/**/*.png
 parser/parser-cp2k/cp2kparser/versions/**/input_data/*.xml
 parser/parser-cp2k/cp2kparser/versions/**/input_data/*.html
 test/unittests/BASIS_SET
diff --git a/parser/parser-cp2k/cp2kparser/versions/cp2k262/commonmatcher.py b/parser/parser-cp2k/cp2kparser/versions/cp2k262/commonmatcher.py
index 0920af7..6a2a4db 100644
--- a/parser/parser-cp2k/cp2kparser/versions/cp2k262/commonmatcher.py
+++ b/parser/parser-cp2k/cp2kparser/versions/cp2k262/commonmatcher.py
@@ -37,31 +37,8 @@ class CommonMatcher(object):
             'x_cp2k_md_force_atom_float': CachingLevel.Cache,
         }
 
-    def adHoc_cp2k_section_cell(self):
-        """Used to extract the cell information.
-        """
-        def wrapper(parser):
-            # Read the lines containing the cell vectors
-            a_line = parser.fIn.readline()
-            b_line = parser.fIn.readline()
-            c_line = parser.fIn.readline()
-
-            # Define the regex that extracts the components and apply it to the lines
-            regex_string = r" CELL\| Vector \w \[angstrom\]:\s+({0})\s+({0})\s+({0})".format(self.regex_f)
-            regex_compiled = re.compile(regex_string)
-            a_result = regex_compiled.match(a_line)
-            b_result = regex_compiled.match(b_line)
-            c_result = regex_compiled.match(c_line)
-
-            # Convert the string results into a 3x3 numpy array
-            cell = np.zeros((3, 3))
-            cell[0, :] = [float(x) for x in a_result.groups()]
-            cell[1, :] = [float(x) for x in b_result.groups()]
-            cell[2, :] = [float(x) for x in c_result.groups()]
-
-            # Push the results to the correct section
-            parser.backend.addArrayValues("simulation_cell", cell, unit="angstrom")
-        return wrapper
+    #===========================================================================
+    # SimpleMatcher trees
 
     # SimpleMatcher for the header that is common to all run types
     def header(self):
@@ -84,6 +61,7 @@ class CommonMatcher(object):
                 ),
                 SM( r" CP2K\| Input file name\s+(?P<x_cp2k_input_filename>.+$)",
                     sections=['x_cp2k_section_filenames'],
+                    otherMetaInfo=['section_XC_functionals', 'XC_functional_name', 'XC_functional_weight', 'XC_functional', 'configuration_periodic_dimensions', "stress_tensor_method"],
                     subMatchers=[
                         SM( r" GLOBAL\| Basis set file name\s+(?P<x_cp2k_basis_set_filename>.+$)"),
                         SM( r" GLOBAL\| Geminal file name\s+(?P<x_cp2k_geminal_filename>.+$)"),
@@ -115,6 +93,8 @@ class CommonMatcher(object):
             ]
         )
 
+    #===========================================================================
+    # Section close triggers
     def onClose_section_method(self, backend, gIndex, section):
         """When all the functional definitions have been gathered, matches them
         with the nomad correspondents and combines into one single string which
@@ -122,19 +102,22 @@ class CommonMatcher(object):
         """
         # Transform the CP2K self-interaction correction string to the NOMAD
         # correspondent, and push directly to the superBackend to avoid caching
-        sic_cp2k = section["self_interaction_correction_method"][0]
-        sic_map = {
-            "NO": "",
-            "AD SIC": "SIC_AD",
-            "Explicit Orbital SIC": "SIC_EXPLICIT_ORBITALS",
-            "SPZ/MAURI SIC": "SIC_MAURI_SPZ",
-            "US/MAURI SIC": "SIC_MAURI_US",
-        }
-        sic_nomad = sic_map.get(sic_cp2k)
-        if sic_nomad is not None:
-            backend.superBackend.addValue('self_interaction_correction_method', sic_nomad)
-        else:
-            logger.warning("Unknown self-interaction correction method used.")
+        try:
+            sic_cp2k = section["self_interaction_correction_method"][0]
+            sic_map = {
+                "NO": "",
+                "AD SIC": "SIC_AD",
+                "Explicit Orbital SIC": "SIC_EXPLICIT_ORBITALS",
+                "SPZ/MAURI SIC": "SIC_MAURI_SPZ",
+                "US/MAURI SIC": "SIC_MAURI_US",
+            }
+            sic_nomad = sic_map.get(sic_cp2k)
+            if sic_nomad is not None:
+                backend.superBackend.addValue('self_interaction_correction_method', sic_nomad)
+            else:
+                logger.warning("Unknown self-interaction correction method used.")
+        except:
+            pass
 
     def onClose_x_cp2k_section_filenames(self, backend, gIndex, section):
         """
@@ -148,6 +131,34 @@ class CommonMatcher(object):
         else:
             logger.warning("The input file of the calculation could not be found.")
 
+    #===========================================================================
+    # Ad hoc parsing
+    def adHoc_cp2k_section_cell(self):
+        """Used to extract the cell information.
+        """
+        def wrapper(parser):
+            # Read the lines containing the cell vectors
+            a_line = parser.fIn.readline()
+            b_line = parser.fIn.readline()
+            c_line = parser.fIn.readline()
+
+            # Define the regex that extracts the components and apply it to the lines
+            regex_string = r" CELL\| Vector \w \[angstrom\]:\s+({0})\s+({0})\s+({0})".format(self.regex_f)
+            regex_compiled = re.compile(regex_string)
+            a_result = regex_compiled.match(a_line)
+            b_result = regex_compiled.match(b_line)
+            c_result = regex_compiled.match(c_line)
+
+            # Convert the string results into a 3x3 numpy array
+            cell = np.zeros((3, 3))
+            cell[0, :] = [float(x) for x in a_result.groups()]
+            cell[1, :] = [float(x) for x in b_result.groups()]
+            cell[2, :] = [float(x) for x in c_result.groups()]
+
+            # Push the results to the correct section
+            parser.backend.addArrayValues("simulation_cell", cell, unit="angstrom")
+        return wrapper
+
     def getOnCloseTriggers(self):
         """
         Returns:
diff --git a/parser/parser-cp2k/cp2kparser/versions/cp2k262/inputparser.py b/parser/parser-cp2k/cp2kparser/versions/cp2k262/inputparser.py
index 6bab307..48dec90 100644
--- a/parser/parser-cp2k/cp2kparser/versions/cp2k262/inputparser.py
+++ b/parser/parser-cp2k/cp2kparser/versions/cp2k262/inputparser.py
@@ -125,6 +125,20 @@ class CP2KInputParser(BasicParser):
             force_file_path = self.normalize_cp2k_path(force_file, "xyz")
             self.file_service.set_file_id(force_file_path, "force_file_single_point")
 
+        #=======================================================================
+        # Stress tensor calculation method
+        stress_tensor_method = self.input_tree.get_keyword("FORCE_EVAL/STRESS_TENSOR")
+        if stress_tensor_method != "NONE":
+            mapping = {
+                "NUMERICAL": "Numerical",
+                "ANALYTICAL": "Analytical",
+                "DIAGONAL_ANALYTICAL": "Diagonal analytical",
+                "DIAGONAL_NUMERICAL": "Diagonal numerical",
+            }
+            stress_tensor_method = mapping.get(stress_tensor_method)
+            if stress_tensor_method is not None:
+                self.backend.addValue("stress_tensor_method", stress_tensor_method)
+
     def normalize_cp2k_path(self, path, extension, name=""):
         """The paths in CP2K input can be given in many ways. This function
         tries to normalize these forms into a valid path.
diff --git a/parser/parser-cp2k/cp2kparser/versions/cp2k262/singlepointparser.py b/parser/parser-cp2k/cp2kparser/versions/cp2k262/singlepointparser.py
index dcdd0bf..0fde035 100644
--- a/parser/parser-cp2k/cp2kparser/versions/cp2k262/singlepointparser.py
+++ b/parser/parser-cp2k/cp2kparser/versions/cp2k262/singlepointparser.py
@@ -43,9 +43,11 @@ class CP2KSinglePointParser(MainHierarchicalParser):
                             ]
                         ),
                         SM( r"  \*\*\* SCF run converged in\s+(\d+) steps \*\*\*",
+                            otherMetaInfo=["single_configuration_calculation_converged"],
                             adHoc=self.adHoc_single_point_converged()
                         ),
                         SM( r"  \*\*\* SCF run NOT converged \*\*\*",
+                            otherMetaInfo=["single_configuration_calculation_converged"],
                             adHoc=self.adHoc_single_point_not_converged()
                         ),
                         SM( r"  Electronic kinetic energy:\s+(?P<electronic_kinetic_energy__hartree>{})".format(self.cm.regex_f)),
@@ -55,12 +57,22 @@ class CP2KSinglePointParser(MainHierarchicalParser):
                         SM( r" ENERGY\| Total FORCE_EVAL \( \w+ \) energy \(a\.u\.\):\s+(?P<energy_total__hartree>{0})".format(self.cm.regex_f)),
                         SM( r" ATOMIC FORCES in \[a\.u\.\]"),
                         SM( r" # Atom   Kind   Element          X              Y              Z",
-                            adHoc=self.adHoc_atom_forces()
-                        ),
-                        SM( r" NUMERICAL STRESS TENSOR [GPa]"),
-                        SM( r"\s+X\s+Y\s+Z",
-                            adHoc=self.adHoc_stress_tensor()
+                            adHoc=self.adHoc_atom_forces(),
+                            otherMetaInfo=["atom_forces"],
                         ),
+                        SM( r" (?:NUMERICAL )?STRESS TENSOR \[GPa\]",
+                            sections=["section_stress_tensor"],
+                            otherMetaInfo=["stress_tensor"],
+                            subMatchers=[
+                                SM( r"\s+X\s+Y\s+Z",
+                                    adHoc=self.adHoc_stress_tensor()
+                                ),
+                                SM( "  1/3 Trace\(stress tensor\):\s+(?P<x_cp2k_stress_tensor_one_third_of_trace__GPa>{})".format(self.cm.regex_f)),
+                                SM( "  Det\(stress tensor\)\s+:\s+(?P<x_cp2k_stress_tensor_determinant__GPa3>{})".format(self.cm.regex_f)),
+                                SM( " EIGENVECTORS AND EIGENVALUES OF THE STRESS TENSOR",
+                                    adHoc=self.adHoc_stress_tensor_eigenpairs()),
+                            ]
+                        )
                     ]
                 )
             ]
@@ -72,6 +84,7 @@ class CP2KSinglePointParser(MainHierarchicalParser):
         self.root_matcher = SM("",
             forwardMatch=True,
             sections=['section_run', "section_single_configuration_calculation", "section_system", "section_method"],
+            otherMetaInfo=["atom_forces"],
             subMatchers=[
                 self.cm.header(),
                 self.energy_force
@@ -238,10 +251,23 @@ class CP2KSinglePointParser(MainHierarchicalParser):
             row2 = [float(x) for x in parser.fIn.readline().split()[-3:]]
             row3 = [float(x) for x in parser.fIn.readline().split()[-3:]]
             stress_array = np.array([row1, row2, row3])
-            gid = parser.backend.openSection("section_stress_tensor")
-            parser.backend.addArrayValues("stress_tensor_value", stress_array, unit="GPa")
-            parser.backend.closeSection("section_stress_tensor", gid)
+            parser.backend.addArrayValues("stress_tensor", stress_array, unit="GPa")
+
+        return wrapper
 
+    def adHoc_stress_tensor_eigenpairs(self):
+        """Parses the stress tensor eigenpairs.
+        """
+        def wrapper(parser):
+            parser.fIn.readline()
+            eigenvalues = np.array([float(x) for x in parser.fIn.readline().split()][::-1])
+            parser.fIn.readline()
+            row1 = [float(x) for x in parser.fIn.readline().split()]
+            row2 = [float(x) for x in parser.fIn.readline().split()]
+            row3 = [float(x) for x in parser.fIn.readline().split()]
+            eigenvectors = np.fliplr(np.array([row1, row2, row3]))
+            parser.backend.addArrayValues("x_cp2k_stress_tensor_eigenvalues", eigenvalues, unit="GPa")
+            parser.backend.addArrayValues("x_cp2k_stress_tensor_eigenvectors", eigenvectors)
         return wrapper
 
     def adHoc_single_point_converged(self):
diff --git a/test/unittests/cp2k_2.6.2/run_tests.py b/test/unittests/cp2k_2.6.2/run_tests.py
index 81f8ec7..9d3c37c 100644
--- a/test/unittests/cp2k_2.6.2/run_tests.py
+++ b/test/unittests/cp2k_2.6.2/run_tests.py
@@ -170,6 +170,39 @@ class TestSelfInteractionCorrectionMethod(unittest.TestCase):
         self.assertEqual(sic, "SIC_MAURI_US")
 
 
+#===============================================================================
+class TestStressTensorMethods(unittest.TestCase):
+    """Tests that the stress tensor can be properly parsed for different
+    calculation methods.
+    """
+    def test_none(self):
+        get_results("stress_tensor/none", "section_stress_tensor")
+
+    def test_analytical(self):
+        results = get_results("stress_tensor/analytical", ["stress_tensor_method", "stress_tensor"])
+        method = results["stress_tensor_method"]
+        results["stress_tensor"]
+        self.assertEqual(method, "Analytical")
+
+    def test_numerical(self):
+        results = get_results("stress_tensor/numerical", ["stress_tensor_method", "stress_tensor"])
+        method = results["stress_tensor_method"]
+        results["stress_tensor"]
+        self.assertEqual(method, "Numerical")
+
+    def test_diagonal_analytical(self):
+        results = get_results("stress_tensor/diagonal_analytical", ["stress_tensor_method", "stress_tensor"])
+        method = results["stress_tensor_method"]
+        results["stress_tensor"]
+        self.assertEqual(method, "Diagonal analytical")
+
+    def test_diagonal_numerical(self):
+        results = get_results("stress_tensor/diagonal_numerical", ["stress_tensor_method", "stress_tensor"])
+        method = results["stress_tensor_method"]
+        results["stress_tensor"]
+        self.assertEqual(method, "Diagonal numerical")
+
+
 #===============================================================================
 class TestConfigurationPeriodicDimensions(unittest.TestCase):
     """Tests that the self-interaction correction can be properly parsed.
@@ -336,7 +369,7 @@ class TestEnergyForce(unittest.TestCase):
         self.assertEqual(result, 0)
 
     def test_stress_tensor(self):
-        result = self.results["stress_tensor_value"]
+        result = self.results["stress_tensor"]
         expected_result = convert_unit(
             np.array([
                 [7.77641684, -0.00000106, -0.00000106],
@@ -347,6 +380,30 @@ class TestEnergyForce(unittest.TestCase):
         )
         self.assertTrue(np.array_equal(result, expected_result))
 
+    def test_stress_tensor_eigenvalues(self):
+        result = self.results["x_cp2k_stress_tensor_eigenvalues"]
+        expected_result = convert_unit(np.array([7.77641809, 7.77641797, 7.77641485]), "GPa")
+        self.assertTrue(np.array_equal(result, expected_result))
+
+    def test_stress_tensor_eigenvectors(self):
+        result = self.results["x_cp2k_stress_tensor_eigenvectors"]
+        expected_result = np.array([
+            [0.00094549, -0.79967815, 0.60042815],
+            [-0.70749682, 0.42379757, 0.56554741],
+            [0.70671590, 0.42533573, 0.56536905],
+        ])
+        self.assertTrue(np.array_equal(result, expected_result))
+
+    def test_stress_tensor_determinant(self):
+        result = self.results["x_cp2k_stress_tensor_determinant"]
+        expected_result = convert_unit(4.70260626E+02, "GPa^3")
+        self.assertTrue(np.array_equal(result, expected_result))
+
+    def test_stress_tensor_one_third_of_trace(self):
+        result = self.results["x_cp2k_stress_tensor_one_third_of_trace"]
+        expected_result = convert_unit(7.77641697E+00, "GPa")
+        self.assertTrue(np.array_equal(result, expected_result))
+
 #===============================================================================
 if __name__ == '__main__':
     pass
@@ -357,6 +414,7 @@ if __name__ == '__main__':
     suites.append(unittest.TestLoader().loadTestsFromTestCase(TestErrors))
     suites.append(unittest.TestLoader().loadTestsFromTestCase(TestXCFunctional))
     suites.append(unittest.TestLoader().loadTestsFromTestCase(TestEnergyForce))
+    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestStressTensorMethods))
     suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSelfInteractionCorrectionMethod))
     suites.append(unittest.TestLoader().loadTestsFromTestCase(TestConfigurationPeriodicDimensions))
     suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSCFConvergence))
-- 
GitLab