adapt abacus/md interface to MD output without stress calculation. (#301)

* adapt abacus/md interface to MD output without stress calculation.

* add unit test for abacus md without stress output.

* Add files via upload

* add assertion failed reason for abacus/md interface

Co-authored-by: LiuRenxi <[email protected]>

Author: Liu-RX

dpdata/abacus/md.py

@@ -33,8 +33,17 @@ def get_coord_dump_freq(inlines):
 def get_coords_from_dump(dumplines, natoms):
     nlines = len(dumplines)
     total_natoms = sum(natoms)
-    nframes_dump = int(nlines/(total_natoms + 13))
-    
+    calc_stress = False
+    if "VIRIAL" in dumplines[6]:
+        calc_stress = True
+    else:
+        assert("POSITIONS" in dumplines[6] and "FORCE" in dumplines[6]), "keywords 'POSITIONS' and 'FORCE' cannot be found in the 6th line. Please check."
+    nframes_dump = -1
+    if calc_stress:
+        nframes_dump = int(nlines/(total_natoms + 13))
+    else:
+        nframes_dump = int(nlines/(total_natoms + 9))
+    assert(nframes_dump > 0), "Number of lines in MD_dump file = %d. Number of atoms = %d. The MD_dump file is incomplete."%(nlines, total_natoms)
     cells = np.zeros([nframes_dump, 3, 3])
     stresses = np.zeros([nframes_dump, 3, 3])
     forces = np.zeros([nframes_dump, total_natoms, 3])
@@ -47,12 +56,17 @@ def get_coords_from_dump(dumplines, natoms):
             # read in LATTICE_VECTORS
             for ix in range(3):
                 cells[iframe, ix] = np.array([float(i) for i in re.split('\s+', dumplines[iline+3+ix])[-3:]]) * celldm
-                stresses[iframe, ix] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+ix])[-3:]])
+                if calc_stress:
+                    stresses[iframe, ix] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+ix])[-3:]])
             for iat in range(total_natoms):
-                coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-6:-3]])*celldm
-                forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-3:]])
+                if calc_stress:
+                    coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-6:-3]])*celldm
+                    forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+11+iat])[-3:]])
+                else:
+                    coords[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+iat])[-6:-3]])*celldm
+                    forces[iframe, iat] = np.array([float(i) for i in re.split('\s+', dumplines[iline+7+iat])[-3:]])
             iframe += 1
-    assert(iframe == nframes_dump)
+    assert(iframe == nframes_dump), "iframe=%d, nframe_dump=%d. Number of frames does not match number of lines in MD_dump."%(iframe, nframes_dump)
     cells *= bohr2ang
     coords *= bohr2ang
     stresses *= kbar2evperang3
@@ -66,7 +80,7 @@ def get_energy(outlines, ndump, dump_freq):
             if nenergy%dump_freq == 0:
                 energy.append(float(line.split()[-2]))
             nenergy+=1
-    assert(ndump == len(energy))
+    assert(ndump == len(energy)), "Number of total energies in running_md.log = %d. Number of frames in MD_dump = %d. Please check."%(len(energy), ndump)
     energy = np.array(energy)
     return energy
 

tests/abacus.md.nostress/INPUT

@@ -0,0 +1,27 @@
+INPUT_PARAMETERS
+#Parameters	(General)
+suffix          autotest
+pseudo_dir		./
+ntype			1	
+nbands			8
+calculation     md
+
+#Parameters (Accuracy)
+ecutwfc			20
+scf_nmax			20
+
+basis_type		pw
+md_nstep           3
+
+#cal_stress          1
+stress_thr      1e-6
+cal_force           1
+force_thr_ev    1.0e-3
+
+ks_solver       cg
+mixing_type     pulay
+mixing_beta     0.7
+
+md_type         -1
+md_dt           1
+md_tfirst       0

tests/abacus.md.nostress/OUT.autotest/MD_dump

@@ -0,0 +1,44 @@
+MDSTEP:  0
+LATTICE_CONSTANT: 10.200000000000
+LATTICE_VECTORS
+  0.500000000000  0.500000000000  0.000000000000
+  0.500000000000  0.000000000000  0.500000000000
+  0.000000000000  0.500000000000  0.500000000000
+INDEX    LABEL    POSITIONS    FORCE (eV/Angstrom)
+  0  Si  0.000000000000  0.000000000000  0.000000000000  -0.885362725259  0.500467424429  0.150239620221
+  1  Si  0.241000000000  0.255000000000  0.250999999999  0.885362725259  -0.500467424429  -0.150239620221
+
+
+MDSTEP:  1
+LATTICE_CONSTANT: 10.200000000000
+LATTICE_VECTORS
+  0.500000000000  0.500000000000  0.000000000000
+  0.500000000000  0.000000000000  0.500000000000
+  0.000000000000  0.500000000000  0.500000000000
+INDEX    LABEL    POSITIONS    FORCE (eV/Angstrom)
+  0  Si  0.999208682339  0.500447306737  0.500134280856  -0.728804779648  0.408578746723  0.107042476463
+  1  Si  0.241791317661  0.254552693263  0.250865719143  0.728804779648  -0.408578746723  -0.107042476463
+
+
+MDSTEP:  2
+LATTICE_CONSTANT: 10.200000000000
+LATTICE_VECTORS
+  0.500000000000  0.500000000000  0.000000000000
+  0.500000000000  0.000000000000  0.500000000000
+  0.000000000000  0.500000000000  0.500000000000
+INDEX    LABEL    POSITIONS    FORCE (eV/Angstrom)
+  0  Si  0.997114226602  0.501624803733  0.500458153134  -0.316038867402  0.171613491764  0.014802659803
+  1  Si  0.243885773398  0.253375196267  0.250541846865  0.316038867402  -0.171613491764  -0.014802659803
+
+
+MDSTEP:  3
+LATTICE_CONSTANT: 10.200000000000
+LATTICE_VECTORS
+  0.500000000000  0.500000000000  0.000000000000
+  0.500000000000  0.000000000000  0.500000000000
+  0.000000000000  0.500000000000  0.500000000000
+INDEX    LABEL    POSITIONS    FORCE (eV/Angstrom)
+  0  Si  0.994451593844  0.503106810892  0.500786060548  0.204324467775  -0.117398116295  -0.052955519257
+  1  Si  0.246548406156  0.251893189108  0.250213939451  -0.204324467775  0.117398116295  0.052955519257
+
+