Bug Fixes: virial should be stress multiplied by volume.

Bug Fixes: virial should be stress multiplied by volume.

cp2kdata/block_parser/cells.py

@@ -95,7 +95,7 @@ def parse_all_md_cells(output_file: List[str],
     # notice that the cell of step 0 is excluded from MD| block
 
     # choose parser according to cp2k_info.version
-    if cp2k_info.version in ['2023.1']:
+    if cp2k_info.version in ['2022.2', '2023.1']:
         ALL_MD_CELL_RE = ALL_MD_CELL_RE_V2023
     elif cp2k_info.version in ['7.1']:
         ALL_MD_CELL_RE = ALL_MD_CELL_RE_V7

cp2kdata/dpdata_plugin.py

@@ -12,6 +12,11 @@
 
 
 WRAPPER = "--- You are parsing data using package Cp2kData ---"
+VIRIAL_WRN = (
+    "Virial Parsing using cp2kdata as plug in for dpdata\n"
+    "was not multiplied by volume before cp2kdata v0.6.4\n"
+    "please check the cp2kdata version and the virial.npy\n"
+    )
 
 
 @Format.register("cp2k/output")
@@ -55,7 +60,10 @@ def from_labeled_system(self, file_name, **kwargs):
         data['coords'] = cp2k_e_f.init_atomic_coordinates[np.newaxis, :, :]
         data['forces'] = cp2k_e_f.atomic_forces_list * AU_TO_EV/AU_TO_ANG
         if cp2k_e_f.has_stress():
-            data['virials'] = cp2k_e_f.stress_tensor_list/EV_ANG_m3_TO_GPa
+            # note that virial = stress * volume
+            print(VIRIAL_WRN)
+            volume = np.linalg.det(data['cells'][0])
+            data['virials'] = cp2k_e_f.stress_tensor_list*volume/EV_ANG_m3_TO_GPa
 
         print(WRAPPER)
         return data
@@ -116,7 +124,13 @@ def from_labeled_system(self, file_name, **kwargs):
         data['coords'] = cp2kmd.atomic_frames_list
         data['forces'] = cp2kmd.atomic_forces_list * AU_TO_EV/AU_TO_ANG
         if cp2kmd.has_stress():
-            data['virials'] = cp2kmd.stress_tensor_list/EV_ANG_m3_TO_GPa
+            # note that virial = stress * volume
+            print(VIRIAL_WRN)
+            # the shape of cells should be (num_frames, 3, 3)
+            # the np.linalg.det() function can handle this and return (num_frames,)
+            volumes = np.linalg.det(data['cells'])
+            volumes = volumes[:, np.newaxis, np.newaxis]
+            data['virials'] = cp2kmd.stress_tensor_list*volumes/EV_ANG_m3_TO_GPa
 
         print(WRAPPER)
         return data

cp2kdata/output.py

@@ -397,7 +397,7 @@ def parse_md(self):
                 "\n"
                 "cp2kdata is parsing md cell information from output file.\n"
                 "The raw data of cell information are lengths and angles,\n"
-                "which are latter transformed to cell matrices by codes.\n"
+                "which are later transformed to cell matrices by codes.\n"
                 "However, the a axis of the cell are always assumed to be aligned to "
                 "the x axis of the coordinate.\n"
                 "Make sure the a axis in real cell matrices are always aligned to x axis.\n"

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "Cp2kData"
-version = "0.6.3"
+version = "0.6.4"
 description = "A Small Package to Postprocess Cp2k Output"
 authors = [
     {name = "Yongbin Zhuang", email = "[email protected]"}

tests/test_dpdata/test_labeledsys.py

@@ -18,6 +18,7 @@
     "tests/test_dpdata/v7.1/aimd_npt_f",
     "tests/test_dpdata/v9.1/xTBmd_npt_i",
     "tests/test_dpdata/v2022.1/aimd",
+    "tests/test_dpdata/v2022.2/aimd_npt_i",
     "tests/test_dpdata/v2023.1/aimd_nvt",
     "tests/test_dpdata/v2023.1/aimd_npt_f"
 ]
@@ -64,7 +65,7 @@
 def cp2k_and_ref(request):
     return request.param
 
-
+#@pytest.mark.parametrize(id=(e_f_output_path_list+aimd_output_path_list))
 class TestLabeledSys():
     def test_len_func(self, cp2k_and_ref):
         assert len(cp2k_and_ref[0]) == len(cp2k_and_ref[1])
@@ -126,9 +127,13 @@ def test_force(self, cp2k_and_ref):
             )
 
     def test_virial(self, cp2k_and_ref):
-        if not 'virials' in cp2k_and_ref[0]:
-            assert False == ('virials' in cp2k_and_ref[1])
+        print(cp2k_and_ref[0])
+        print(cp2k_and_ref[0].has_virial())
+        if not cp2k_and_ref[0].has_virial():
+
+            assert False == (cp2k_and_ref[1].has_virial())
             return
+        print("if code works ok")
         np.testing.assert_almost_equal(
             cp2k_and_ref[0]['virials'],
             cp2k_and_ref[1]['virials'],

tests/test_dpdata/v2022.2/aimd_npt_i/coord.xyz

@@ -0,0 +1,52 @@
+Li  2.46250000  2.46250000  0.22766305
+Li  2.46250000  2.46250000  4.69733695
+Li  2.46250000  4.69733695  2.46250000
+Li  2.46250000  0.22766305  2.46250000
+Li  0.22766305  2.46250000  2.46250000
+Li  4.69733695  2.46250000  2.46250000
+Li  2.46250000  7.38750000  5.15266305
+Li  2.46250000  7.38750000  9.62233695
+Li  2.46250000  9.62233695  7.38750000
+Li  2.46250000  5.15266305  7.38750000
+Li  0.22766305  7.38750000  7.38750000
+Li  4.69733695  7.38750000  7.38750000
+Li  7.38750000  2.46250000  5.15266305
+Li  7.38750000  2.46250000  9.62233695
+Li  7.38750000  4.69733695  7.38750000
+Li  7.38750000  0.22766305  7.38750000
+Li  5.15266305  2.46250000  7.38750000
+Li  9.62233695  2.46250000  7.38750000
+Li  7.38750000  7.38750000  0.22766305
+Li  7.38750000  7.38750000  4.69733695
+Li  7.38750000  9.62233695  2.46250000
+Li  7.38750000  5.15266305  2.46250000
+Li  5.15266305  7.38750000  2.46250000
+Li  9.62233695  7.38750000  2.46250000
+P  4.92500000  0.00000000  0.00000000
+P  4.92500000  4.92500000  4.92500000
+P  0.00000000  0.00000000  4.92500000
+P  0.00000000  4.92500000  0.00000000
+S  1.13789495  3.78710505  8.71210505
+S  8.71210505  3.78710505  1.13789495
+S  3.78710505  1.13789495  8.71210505
+S  6.06289495  1.13789495  1.13789495
+S  1.13789495  8.71210505  3.78710505
+S  8.71210505  8.71210505  6.06289495
+S  3.78710505  6.06289495  3.78710505
+S  6.06289495  6.06289495  6.06289495
+S  6.06289495  3.78710505  3.78710505
+S  3.78710505  3.78710505  6.06289495
+S  8.71210505  1.13789495  3.78710505
+S  1.13789495  1.13789495  6.06289495
+S  6.06289495  8.71210505  8.71210505
+S  3.78710505  8.71210505  1.13789495
+S  8.71210505  6.06289495  8.71210505
+S  1.13789495  6.06289495  1.13789495
+S  7.38750000  2.46250000  7.38750000
+S  7.38750000  7.38750000  2.46250000
+S  0.00000000  0.00000000  0.00000000
+S  0.00000000  4.92500000  4.92500000
+Cl  4.92500000  0.00000000  4.92500000
+Cl  4.92500000  4.92500000  0.00000000
+Cl  2.46250000  2.46250000  2.46250000
+Cl  2.46250000  7.38750000  7.38750000

tests/test_dpdata/v2022.2/aimd_npt_i/deepmd/type.raw

@@ -0,0 +1,52 @@
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+0
+1
+1
+1
+1
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+2
+3
+3
+3
+3

tests/test_dpdata/v2022.2/aimd_npt_i/deepmd/type_map.raw

@@ -0,0 +1,4 @@
+Li
+P
+S
+Cl

tests/test_dpdata/v2022.2/aimd_npt_i/input.inp

@@ -0,0 +1,114 @@
+&GLOBAL
+	PROJECT system
+	RUN_TYPE MD
+	PRINT_LEVEL MEDIUM
+&END
+
+&MOTION
+
+  &MD
+    ENSEMBLE NPT_I
+    STEPS 10
+    TIMESTEP 2.0
+    TEMPERATURE 300
+
+    &BAROSTAT
+      PRESSURE 1
+      TIMECON  100
+    &END BAROSTAT
+
+    &THERMOSTAT
+      &NOSE
+        LENGTH 3
+        YOSHIDA 3
+        TIMECON 200.0
+        MTS 2
+      &END NOSE
+    &END THERMOSTAT
+  &END MD
+
+&END
+
+
+&FORCE_EVAL
+	METHOD Quickstep
+	STRESS_TENSOR ANALYTICAL
+	&PRINT
+		&FORCES ON
+		&END
+		&STRESS_TENSOR ON
+                &END
+	&END
+
+	&DFT
+		BASIS_SET_FILE_NAME  BASIS_MOLOPT
+		POTENTIAL_FILE_NAME  GTH_POTENTIALS
+#		MULTIPLICITY 2
+		&QS
+			EPS_DEFAULT 1.0E-12
+		&END
+		&MGRID
+			NGRIDS 4
+			CUTOFF 1000
+			REL_CUTOFF 60
+		&END
+		&XC
+			&XC_FUNCTIONAL PBE
+			&END
+		&END
+
+		&SCF
+			SCF_GUESS ATOMIC
+			EPS_SCF 1.0E-6
+			MAX_SCF 100
+			ADDED_MOS 200
+			&SMEAR ON
+				METHOD FERMI_DIRAC
+				#ELECTRONIC_TEMPERATURE [K] 300
+			&END SMEAR
+			&DIAGONALIZATION  ON
+				ALGORITHM STANDARD
+			&END 
+			&MIXING  T
+				METHOD BROYDEN_MIXING
+				ALPHA 0.4
+				NBROYDEN 8
+			&END
+		&END 
+		&KPOINTS
+				SCHEME  MONKHORST-PACK 2 2 2
+		&END
+	&END
+
+	&SUBSYS
+		&KIND Li
+			ELEMENT   Li
+			BASIS_SET DZVP-MOLOPT-SR-GTH
+			POTENTIAL GTH-PBE
+		&END
+		&KIND P
+			ELEMENT   P
+			BASIS_SET DZVP-MOLOPT-SR-GTH
+			POTENTIAL GTH-PBE
+		&END
+		&KIND S
+			ELEMENT   S
+			BASIS_SET DZVP-MOLOPT-SR-GTH
+			POTENTIAL GTH-PBE
+		&END
+		&KIND Cl
+			ELEMENT   Cl
+			BASIS_SET DZVP-MOLOPT-SR-GTH
+			POTENTIAL GTH-PBE
+		&END
+		&CELL
+			A 9.8500000000 0.0000000000 0.0000000000 
+			B 0.0000000000 9.8500000000 0.0000000000 
+			C 0.0000000000 0.0000000000 9.8500000000 
+		&END
+
+		&COORD
+		@include coord.xyz
+		&END
+	&END
+&END