Fix virial

dpdata/cp2k/output.py

@@ -120,7 +120,22 @@ def handle_single_log_frame(self, lines):
         print_level_flag = 0
         atomic_kinds_pattern = re.compile(r'\s+\d+\. Atomic kind:\s+(?P<akind>\S+)')
         atomic_kinds = [] 
+        stress_sign = 'STRESS'
+        stress_flag = 0
+        stress = []
+
         for line in lines:
+            if stress_flag == 3 :
+                if (line == '\n') :
+                    stress_flag = 0
+                else :
+                    stress.append(line.split()[1:4])
+            if stress_flag == 2 :
+                stress_flag = 3
+            if stress_flag == 1 :
+                stress_flag = 2
+            if (stress_sign in line):
+                stress_flag = 1
             if force_start_pattern.match(line):
                 force_flag=True
             if force_end_pattern.match(line):
@@ -214,7 +229,18 @@ def handle_single_log_frame(self, lines):
         #atom_names=list(element_dict.keys())
         atom_names=self.atomic_kinds
         atom_numbs=[]
-
+
+        GPa = PressureConversion("eV/angstrom^3", "GPa").value()
+        if stress:
+            stress = np.array(stress)
+            stress = stress.astype('float32')
+            stress = stress[np.newaxis, :, :]
+            # stress to virial conversion, default unit in cp2k is GPa
+            # note the stress is virial = stress * volume
+            virial = stress * np.linalg.det(self.cell)/GPa
+            virial = virial.squeeze()
+        else:
+            virial = None
         for ii in element_dict.keys():
             atom_numbs.append(element_dict[ii][1])
         #print(atom_numbs)
@@ -225,6 +251,8 @@ def handle_single_log_frame(self, lines):
         info_dict['cells'] = np.asarray([self.cell]).astype('float32')
         info_dict['energies'] = np.asarray([energy]).astype('float32')
         info_dict['forces'] = np.asarray([forces_list]).astype('float32')
+        if(virial is not None ): 
+            info_dict['virials'] = np.asarray([virial]).astype('float32')
         return info_dict
 
     def handle_single_xyz_frame(self, lines):

tests/cp2k/aimd_stress/DPGEN-pos-1.xyz

@@ -0,0 +1,7 @@
+     5
+ i =        0, time =        0.000, E =      -8.07218972206
+  H 5.70191 3.8096  4.38845
+  H 4.44601 4.65681 5.38195
+  H 4.49216 4.90991 3.57756
+  H 5.83825 5.62228 4.61068
+  C 5.11666 4.75533 4.46123

tests/cp2k/aimd_stress/deepmd/box.raw

@@ -0,0 +1 @@
+1.000000000000000000e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.000000000000000000e+01 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 1.000000000000000000e+01

tests/cp2k/aimd_stress/deepmd/coord.raw

@@ -0,0 +1 @@
+5.701910018920898438e+00 3.809600114822387695e+00 4.388450145721435547e+00 4.446010112762451172e+00 4.656809806823730469e+00 5.381949901580810547e+00 4.492159843444824219e+00 4.909910202026367188e+00 3.577559947967529297e+00 5.838250160217285156e+00 5.622280120849609375e+00 4.610680103302001953e+00 5.116660118103027344e+00 4.755330085754394531e+00 4.461229801177978516e+00

tests/cp2k/aimd_stress/deepmd/energy.raw

@@ -0,0 +1 @@
+-2.196554718017578125e+02

tests/cp2k/aimd_stress/deepmd/force.raw

@@ -0,0 +1 @@
+-2.743706703186035156e-01 4.712709188461303711e-01 -1.367179006338119507e-01 7.652513980865478516e-01 2.514204978942871094e-01 -1.059429287910461426e+00 2.538656070828437805e-02 1.183478906750679016e-02 -1.267343163490295410e-01 -8.532900810241699219e-01 -9.029597043991088867e-01 -1.871660351753234863e-01 3.367811143398284912e-01 1.687173396348953247e-01 1.509592533111572266e+00

tests/cp2k/aimd_stress/deepmd/type.raw

@@ -0,0 +1,5 @@
+0
+0
+0
+0
+1

tests/cp2k/aimd_stress/deepmd/type_map.raw

@@ -0,0 +1,2 @@
+H
+C

tests/cp2k/aimd_stress/deepmd/virial.raw

@@ -0,0 +1 @@
+-1.305763244628906250e+00 -5.517681837081909180e-01 5.745944380760192871e-01 -5.517681837081909180e-01 -1.251762151718139648e+00 5.180240049958229065e-02 5.745944380760192871e-01 5.180240049958229065e-02 -8.817730545997619629e-01

tests/cp2k/aimd_stress/input.inp

@@ -0,0 +1,108 @@
+
+&GLOBAL
+  PRINT_LEVEL medium
+  WALLTIME 43000
+  PROJECT_NAME DPGEN
+  RUN_TYPE MD
+&END GLOBAL
+&MOTION
+  &MD
+    &THERMOSTAT
+      TYPE NOSE
+    &END THERMOSTAT
+    ENSEMBLE NVT
+    STEPS 3
+    TIMESTEP 0.5
+    TEMPERATURE 330.0
+  &END MD
+  &PRINT
+    &FORCES
+      &EACH
+        MD 1
+      &END EACH
+    &END FORCES
+    &TRAJECTORY
+      &EACH
+        MD 1
+      &END EACH
+    &END TRAJECTORY
+    &RESTART
+      BACKUP_COPIES 3
+      &EACH
+        MD 1
+      &END EACH
+    &END RESTART
+  &END PRINT
+&END MOTION
+&FORCE_EVAL
+  METHOD QS
+  STRESS_TENSOR ANALYTICAL
+  &PRINT
+    &FORCES ON
+    &END FORCES
+    &STRESS_TENSOR ON
+    &END STRESS_TENSOR
+  &END PRINT
+  &SUBSYS
+    &CELL
+      A 10.  0.  0.
+      B  0. 10.  0.
+      C  0.  0. 10.
+    &END CELL
+    &COORD
+      H 5.70191 3.8096  4.38845
+      H 4.44601 4.65681 5.38195
+      H 4.49216 4.90991 3.57756
+      H 5.83825 5.62228 4.61068
+      C 5.11666 4.75533 4.46123
+    &END COORD
+    &KIND H
+      POTENTIAL GTH-PBE-q1
+      BASIS_SET DZVP-GTH-PBE
+    &END KIND
+    &KIND C
+      POTENTIAL GTH-PBE-q4
+      BASIS_SET DZVP-GTH-PBE
+    &END KIND
+  &END SUBSYS
+  &DFT
+    POTENTIAL_FILE_NAME GTH_POTENTIALS
+    PLUS_U_METHOD MULLIKEN
+    UKS .TRUE.
+    BASIS_SET_FILE_NAME BASIS_SET
+    &MGRID
+      REL_CUTOFF 60
+      NGRIDS 5
+      CUTOFF 800
+    &END MGRID
+    &QS
+      EPS_DEFAULT 1.0E-13
+    &END QS
+    &XC
+      &VDW_POTENTIAL
+        POTENTIAL_TYPE PAIR_POTENTIAL
+        &PAIR_POTENTIAL
+          REFERENCE_FUNCTIONAL PBE
+          PARAMETER_FILE_NAME dftd3.dat
+          TYPE DFTD3
+        &END PAIR_POTENTIAL
+      &END VDW_POTENTIAL
+      &XC_FUNCTIONAL PBE
+      &END XC_FUNCTIONAL
+    &END XC
+    &SCF
+      EPS_SCF 1e-06
+      MAX_SCF 50
+      SCF_GUESS RESTART
+      &OUTER_SCF
+        EPS_SCF 1e-06
+        MAX_SCF 15
+      &END OUTER_SCF
+      &OT
+        ENERGY_GAP 0.1
+        PRECONDITIONER FULL_SINGLE_INVERSE
+        MINIMIZER CG
+      &END OT
+    &END SCF
+  &END DFT
+&END FORCE_EVAL

tests/test_cp2k_aimd_output.py

@@ -15,6 +15,15 @@ def setUp(self):
         self.f_places = 6
         self.v_places = 4
 
+class TestCp2kAimdStressOutput(unittest.TestCase, CompLabeledSys):
+    def setUp(self):
+        self.system_1 = dpdata.LabeledSystem('cp2k/aimd_stress',fmt='cp2k/aimd_output')
+        self.system_2 = dpdata.LabeledSystem('cp2k/aimd_stress/deepmd', fmt='deepmd/raw')
+        self.places = 6
+        self.e_places = 6
+        self.f_places = 6
+        self.v_places = 4
+
 #class TestCp2kAimdRestartOutput(unittest.TestCase, CompLabeledSys):
 #    def setUp(self):
 #        self.system_1 = dpdata.LabeledSystem('cp2k/restart_aimd',fmt='cp2k/aimd_output', restart=True)