Merge pull request #51 from robinzyb/devel

parse vibrational frequencies from vibrational analysis output

cp2kdata/block_parser/dft_plus_u.py

@@ -2,7 +2,7 @@
 import numpy as np
 
 PLUS_U_RE = re.compile(
-    """
+    r"""
     \sDFT\+U\soccupations\sof\sspin\s(?P<spin>\d)
     \sfor\sthe\satoms\sof\satomic\skind\s(?P<kind>\d+):\s\w+\s*\n
     \n

cp2kdata/block_parser/hirshfeld.py

@@ -2,7 +2,7 @@
 import numpy as np
 
 HIRSHFELD_RE = re.compile(
-    """
+    r"""
     \s+Hirshfeld\sCharges\s*\n
     \n
     \s\#.+\n

cp2kdata/block_parser/vibration.py

@@ -0,0 +1,23 @@
+import regex as re
+import numpy as np
+
+VIB_FREQ_RE = re.compile(
+    r"""
+    \sVIB\|Frequency\s\(cm\^-1\)(\s{0,15}(?P<vib_freq>[\s-]\d+\.\d+)){1,3}
+    """,
+    re.VERBOSE
+)
+
+
+def parse_vibration_freq_list(output_file):
+
+    vib_freq_list = []
+    for match in VIB_FREQ_RE.finditer(output_file):
+        #vib_freq_list.append(match["energy"])
+        for vib_freq in match.captures('vib_freq'):
+            vib_freq_list.append(vib_freq)
+
+    if vib_freq_list:
+        return np.array(vib_freq_list, dtype=float)
+    else:
+        return None

cp2kdata/output.py

@@ -21,6 +21,7 @@
 from cp2kdata.block_parser.stress import parse_stress_tensor_list
 from cp2kdata.block_parser.cells import parse_all_cells, parse_all_md_cells
 from cp2kdata.block_parser.md_xyz import parse_md_ener, parse_pos_xyz, parse_frc_xyz, parse_md_stress, parse_md_cell
+from cp2kdata.block_parser.vibration import parse_vibration_freq_list
 
 
 class Cp2kOutput:
@@ -104,12 +105,16 @@ def __init__(
             "ENERGY_FORCE": self.parse_energy_force,
             "GEO_OPT": self.parse_geo_opt,
             "CELL_OPT": self.parse_cell_opt,
-            "MD": self.parse_md
+            "MD": self.parse_md,
+            "VIBRATIONAL_ANALYSIS": self.parse_vibrational_analysis
         }
 
         # call corresponding parser for run types
-        parse_run_type = run_type_parser_candidates[self.global_info.run_type]
-        parse_run_type()
+        parse_run_type = run_type_parser_candidates.get(self.global_info.run_type, None)
+        if parse_run_type:
+            parse_run_type()
+        else:
+            f"parser for run type {self.global_info.run_type} is not implemented yet!"
 
         # self.errors_info = parse_errors(self.output_file)
         # if ignore_error:
@@ -337,6 +342,9 @@ def parse_energy_force(self):
         self.atomic_forces_list = parse_atomic_forces_list(self.output_file)
         self.stress_tensor_list = parse_stress_tensor_list(self.output_file)
 
+    def parse_vibrational_analysis(self):
+        self.parse_energy_force()
+
     def parse_geo_opt(self):
         self.init_atomic_coordinates, self.atom_kind_list, self.chemical_symbols = parse_init_atomic_coordinates(
             self.output_file)
@@ -506,6 +514,15 @@ def organize_md_cell(self):
             self.all_cells = np.insert(
                 self.all_cells, 0, first_cell[0], axis=0)
 
+    @cached_property
+    def vib_freq_list(self):
+        assert self.global_info.run_type == "VIBRATIONAL_ANALYSIS", "vibrational frequency is only available for VIBRATIONAL_ANALYSIS run type."
+        # use cached property to prase only once
+        return parse_vibration_freq_list(self.output_file)
+
+    def get_vib_freq_list(self):
+        return self.vib_freq_list
+
     @staticmethod
     def drop_last_info(cp2k_info, array):
         # drop last info parsed from output if it is terminated by request (touch EXIT)
@@ -527,7 +544,7 @@ def get_global_info(run_type=None, filename=None):
     @staticmethod
     def check_run_type(run_type):
         implemented_run_type_parsers = \
-            ["ENERGY_FORCE", "ENERGY", "MD", "GEO_OPT", "CELL_OPT"]
+            ["ENERGY_FORCE", "ENERGY", "MD", "GEO_OPT", "CELL_OPT", "VIBRATIONAL_ANALYSIS"]
         if run_type not in implemented_run_type_parsers:
             raise ValueError(
                 f"Parser for Run Type {run_type} haven't been implemented yet!"

docs/output.md

@@ -92,4 +92,75 @@ Alternatively, you may parse without standard outputs. Consequently, you will lo
 from cp2kdata import Cp2kOutput
 cp2k_output_file = "output_md"
 cp2koutput=Cp2kOutput(run_type="md")
-```
+```
+
+
+## Parse VIBRATIONAL_ANALYSIS outputs
+
+
+```python
+from cp2kdata import Cp2kOutput
+cp2k_output_file = "output"
+cp2koutput=Cp2kOutput(cp2k_output_file)
+```
+```
+
+Cp2k Output Summary
+
+--------------------------------------
+
+Cp2k Version       : 7.1
+
+Run Type           : VIBRATIONAL_ANALYSIS
+
+Atom Numbers       : 84
+
+Frame Numbers      : 1
+
+Force in Output    : No
+
+Stress in Output   : No
+
+Element List       : Ir   O    H
+
+Element Numb       : 20   48   16
+--------------------------------------
+```
+
+The frequencies of normal modes can be easily obtained through
+```python
+
+vib_freq = cp2koutput.get_vib_freq_list()
+print(vib_freq)
+```
+
+```python
+
+array([  86.376632,   97.513035,  121.216151,  129.190651,  132.227782,
+        136.462792,  142.615702,  147.191117,  159.80634 ,  162.154751,
+        168.906485,  170.829935,  173.900984,  177.693402,  183.401291,
+        184.730182,  187.076345,  189.403066,  191.336104,  192.041037,
+        193.255332,  195.888228,  198.675828,  199.922815,  206.417016,
+        210.953256,  218.001682,  219.425306,  222.52259 ,  224.924497,
+        234.303527,  234.913409,  238.023057,  248.354352,  255.550832,
+        258.840498,  261.863363,  276.005055,  283.689727,  291.995009,
+        301.789878,  310.297048,  323.748766,  329.656908,  346.17215 ,
+        353.855895,  360.270334,  376.393988,  378.928125,  392.873287,
+        396.863859,  402.392675,  406.699951,  413.228256,  416.512902,
+        427.329694,  435.765748,  438.283727,  452.313934,  458.465777,
+        460.79924 ,  465.321668,  472.298302,  475.940439,  483.095673,
+        484.024696,  485.946038,  488.381279,  501.80579 ,  508.714296,
+        515.24576 ,  524.550782,  548.347144,  553.286646,  561.806014,
+        568.374944,  575.581685,  581.82834 ,  593.158413,  603.508201,
+        612.45818 ,  615.321528,  617.676424,  628.791272,  632.70203 ,
+        634.782343,  637.505407,  646.571283,  647.068464,  665.866031,
+        668.360042,  674.062077,  675.781939,  691.829219,  692.289219,
+        702.477412,  703.782994,  709.119554,  717.545042,  725.023678,
+        725.974217,  732.039541,  735.313457,  738.154712,  741.528806,
+        743.384134,  752.69417 ,  753.024109,  754.98813 ,  758.675974,
+        760.589761,  767.204112,  773.037322,  778.035202,  790.192429,
+        793.116029,  793.676956,  797.698711,  799.368024,  803.600176,
+        809.272806,  822.118902,  840.188822,  850.058458,  895.024762,
+        905.977585,  973.556033,  992.515143, 1014.623025, 1017.337562,
+       1023.263474, 1027.645573])
+```

tests/test_vibrational_analysis/test_vibrational_analysis.py

@@ -0,0 +1,61 @@
+from cp2kdata.output import Cp2kOutput
+import os
+import json
+import pytest
+import numpy as np
+
+
+def json_to_dict(json_file):
+    with open(json_file, "r") as fp:
+        dict_content = json.load(fp)
+    return dict_content
+
+
+vibrational_analysis_output_path_list = [
+    'tests/test_vibrational_analysis/v7.1/normal',
+
+]
+
+vibrational_analysis_output_list = [
+    Cp2kOutput(os.path.join(path, 'output')) for path in vibrational_analysis_output_path_list
+]
+
+
+test_params = list(zip(vibrational_analysis_output_list,
+                   vibrational_analysis_output_path_list))
+
+
+@pytest.fixture(params=test_params, scope='class', ids=vibrational_analysis_output_path_list)
+def output_and_answer_path(request):
+    return request.param
+
+# @pytest.fixture
+# def answer_in_dict(output_and_answer_path):
+#     answer = json_to_dict(
+#             os.path.join(
+#                 output_and_answer_path[1],
+#                 "answer",
+#                 "answer.json"
+#             )
+#         )
+#     return answer
+
+class TestEnergyForce():
+    def test_run_type(self, output_and_answer_path):
+        a_test_output = output_and_answer_path[0]
+        run_type = a_test_output.get_run_type()
+        assert run_type == 'VIBRATIONAL_ANALYSIS'
+
+    def test_atomic_forces_list(self, output_and_answer_path):
+        a_test_output = output_and_answer_path[0]
+        vib_freq = a_test_output.get_vib_freq_list()
+        if a_test_output.has_force():
+            answer = np.load(
+                os.path.join(
+                    output_and_answer_path[1],
+                    "answer",
+                    "vib_freq.npy"
+                )
+            )
+            assert np.all(vib_freq == answer)
+