adapt dpdata interface to ABACUS 2.2.0, fixed issue#268 about qe/pw/scf

dpdata/abacus/md.py

@@ -1,6 +1,8 @@
+from ast import dump
 import os,sys
 import numpy as np
-from .scf import ry2ev, kbar2evperang3, get_block, get_geometry_in, get_cell, get_coords
+from .scf import ry2ev, bohr2ang, kbar2evperang3, get_block, get_geometry_in, get_cell, get_coords
+import re
 
 # Read in geometries from an ABACUS MD trajectory.
 # The atomic coordinates are read in from generated files in OUT.XXXX.
@@ -24,126 +26,49 @@ def get_path_out(fname, inlines):
 
 def get_coord_dump_freq(inlines):
     for line in inlines:
-        if  len(line)>0 and "md_dumpmdfred" in line and "md_dumpmdfred" == line.split()[0]:
+        if  len(line)>0 and "md_dumpfreq" in line and "md_dumpfreq" == line.split()[0]:
             return int(line.split()[1])
     return 1
 
-# set up a cell according to cell info in cif file.
-# maybe useful later
-'''
-def setup_cell(a, b, c, alpha, beta, gamma):
-    cell = np.zeros(3, 3)
-    cell[0, 0] = a
-    cell[1, 0] = b*np.cos(gamma/180*np.pi)
-    cell[1, 1] = b*np.sin(gamma/180*np.pi)
-    cell[2, 0] = c*np.cos(beta/180*np.pi)
-    cell[2, 1] = c*(b*np.cos(alpha/180*np.pi) - cell[1, 0]*np.cos(beta/180*np.pi))/cell[1, 1]
-    cell[2, 2] = np.sqrt(c**2 - cell[2, 0]**2 - cell[2, 1]**2)
-    return cell
-'''
-
-def get_single_coord_from_cif(pos_file, atom_names, natoms, cell):
-    assert(len(atom_names) == len(natoms))
-    nele = len(atom_names)
+def get_coords_from_dump(dumplines, natoms):
+    nlines = len(dumplines)
     total_natoms = sum(natoms)
-    coord = np.zeros([total_natoms, 3])
-    a = 0
-    b = 0
-    c = 0
-    alpha = 0
-    beta = 0
-    gamma = 0
-    with open(pos_file, "r") as fp:
-        lines = fp.read().split("\n")
-    for line in lines:
-        if "_cell_length_a" in line:
-            a = float(line.split()[1])
-        if "_cell_length_b" in line:
-            b = float(line.split()[1])
-        if "_cell_length_c" in line:
-            c = float(line.split()[1])  
-        if "_cell_angle_alpha" in line:
-            alpha = float(line.split()[1])
-        if "_cell_angle_beta" in line:
-            beta = float(line.split()[1])
-        if "_cell_angle_gamma" in line:
-            gamma = float(line.split()[1])
-    assert(a > 0 and b > 0 and c > 0 and alpha > 0 and beta > 0 and gamma > 0)
-    #cell = setup_cell(a, b, c, alpha, beta, gamma)
-    coord_lines = get_block(lines=lines, keyword="_atom_site_fract_z", skip=0, nlines = total_natoms)
-
-    ia_idx = 0
-    for it in range(nele):
-        for ia in range(natoms[it]):
-            coord_line = coord_lines[ia_idx].split()
-            assert(coord_line[0] == atom_names[it])
-            coord[ia_idx, 0] = float(coord_line[1])
-            coord[ia_idx, 1] = float(coord_line[2])
-            coord[ia_idx, 2] = float(coord_line[3])
-            ia_idx+=1
-    coord = np.matmul(coord, cell)
-    # important! Coordinates are converted to Cartesian coordinate.
-    return coord
+    nframes_dump = int(nlines/(total_natoms + 13))
 
-
-def get_coords_from_cif(ndump, dump_freq, atom_names, natoms, types, path_out, cell):
-    total_natoms = sum(natoms)
-    #cell = np.zeros(ndump, 3, 3)
-    coords = np.zeros([ndump, total_natoms, 3])
-    pos_file = os.path.join(path_out, "STRU_READIN_ADJUST.cif")
-    # frame 0 file is different from any other frames
-    coords[0] = get_single_coord_from_cif(pos_file, atom_names, natoms, cell)
-    for dump_idx in range(1, ndump):
-        pos_file = os.path.join(path_out, "md_pos_%d.cif" %(dump_idx*dump_freq))
-        #print("dump_idx = %s" %dump_idx)
-        coords[dump_idx] = get_single_coord_from_cif(pos_file, atom_names, natoms, cell)
-    return coords
+    cells = np.zeros([nframes_dump, 3, 3])
+    stresses = np.zeros([nframes_dump, 3, 3])
+    forces = np.zeros([nframes_dump, total_natoms, 3])
+    coords = np.zeros([nframes_dump, total_natoms, 3])
+    iframe = 0
+    for iline in range(nlines):
+        if "MDSTEP" in dumplines[iline]:
+            # read in LATTICE_CONSTANT
+            celldm = float(dumplines[iline+1].split(" ")[-1])
+            # 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:]])
+            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:]])
+            iframe += 1
+    assert(iframe == nframes_dump)
+    cells *= bohr2ang
+    coords *= bohr2ang
+    stresses *= kbar2evperang3
+    return coords, cells, forces, stresses
 
-def get_energy_force_stress(outlines, inlines, dump_freq, ndump, natoms, atom_names):
-    stress = None
-    total_natoms = sum(natoms)
-    for line in inlines:
-        if len(line)>0 and "stress" in line and "stress" == line.split()[0] and "1" == line.split()[1]:
-            stress = np.zeros([ndump, 3, 3])
-            break
-    if type(stress) != np.ndarray:
-        print("The ABACUS program has no stress output. Stress will not be read.")
+def get_energy(outlines, ndump, dump_freq):
+    energy = []
     nenergy = 0
-    nforce = 0
-    nstress = 0
-    energy = np.zeros(ndump)
-    force = np.zeros([ndump, total_natoms, 3])
-
     for line_idx, line in enumerate(outlines):
         if "final etot is" in line:
             if nenergy%dump_freq == 0:
-                energy[int(nenergy/dump_freq)] = float(line.split()[-2])
+                energy.append(float(line.split()[-2]))
             nenergy+=1
-        if "TOTAL-FORCE (eV/Angstrom)" in line:
-            for iatom in range(0, total_natoms):
-                force_line = outlines[line_idx+5+iatom]
-                atom_force = [float(i) for i in force_line.split()[1:]]
-                assert(len(atom_force) == 3)
-                atom_force = np.array(atom_force)
-                if nforce%dump_freq == 0:
-                    force[int(nforce/dump_freq), iatom] = atom_force
-            nforce+=1
-            assert(nforce==nenergy)
-        if "TOTAL-STRESS (KBAR)" in line:
-            for idx in range(0, 3):
-                stress_line = outlines[line_idx+4+idx]
-                single_stress = [float(i) for i in stress_line.split()]
-                if len(single_stress) != 3:
-                    print(single_stress)
-                assert(len(single_stress) == 3)
-                single_stress = np.array(single_stress)
-                if nstress%dump_freq == 0:
-                    stress[int(nstress/dump_freq), idx] = single_stress
-            nstress+=1
-            assert(nstress==nforce)
-    if type(stress) == np.ndarray:
-        stress *= kbar2evperang3
-    return energy, force, stress
+    assert(ndump == len(energy))
+    energy = np.array(energy)
+    return energy
 
 
 def get_frame (fname):
@@ -164,23 +89,27 @@ def get_frame (fname):
     atom_names, natoms, types, coords = get_coords(celldm, cell, geometry_inlines, inlines) 
     # This coords is not to be used.
     dump_freq = get_coord_dump_freq(inlines = inlines)
-    ndump = int(os.popen("ls -l %s | grep 'md_pos_' | wc -l" %path_out).readlines()[0])
+    #ndump = int(os.popen("ls -l %s | grep 'md_pos_' | wc -l" %path_out).readlines()[0])
     # number of dumped geometry files
-    coords = get_coords_from_cif(ndump, dump_freq, atom_names, natoms, types, path_out, cell)
-
-    # TODO: Read in energies, forces and pressures.
+    #coords = get_coords_from_cif(ndump, dump_freq, atom_names, natoms, types, path_out, cell)
+    with open(os.path.join(path_out, "MD_dump"), 'r') as fp:
+        dumplines = fp.read().split('\n')
+    coords, cells, force, stress = get_coords_from_dump(dumplines, natoms)
+    ndump = np.shape(coords)[0]
     with open(os.path.join(path_out, "running_md.log"), 'r') as fp:
         outlines = fp.read().split('\n')
-    energy, force, stress = get_energy_force_stress(outlines, inlines, dump_freq, ndump, natoms, atom_names)
-    if type(stress) == np.ndarray:
-        stress *= np.linalg.det(cell)
+    energy = get_energy(outlines, ndump, dump_freq)
+    for iframe in range(ndump):
+        stress[iframe] *= np.linalg.det(cells[iframe, :, :].reshape([3, 3]))
+    if np.sum(np.abs(stress[0])) < 1e-10:
+        stress = None
     data = {}
     data['atom_names'] = atom_names
     data['atom_numbs'] = natoms
     data['atom_types'] = types
-    data['cells'] = np.zeros([ndump, 3, 3])
-    for idx in range(ndump):
-        data['cells'][:, :, :] = cell
+    data['cells'] = cells
+    #for idx in range(ndump):
+    #    data['cells'][:, :, :] = cell
     data['coords'] = coords
     data['energies'] = energy
     data['forces'] = force

dpdata/abacus/scf.py

@@ -16,15 +16,18 @@ def get_block (lines, keyword, skip = 0, nlines = None):
             found = True
             blk_idx = idx + 1 + skip
             line_idx = 0
-            while len(lines[blk_idx]) == 0:
+            while len(lines[blk_idx].split("\s+")) == 0:
                 blk_idx += 1
-            while len(lines[blk_idx]) != 0 and line_idx < nlines and blk_idx != len(lines):
+            while line_idx < nlines and blk_idx != len(lines):
+                if len(lines[blk_idx].split("\s+")) == 0 or lines[blk_idx] == "":
+                    blk_idx+=1
+                    continue
                 ret.append(lines[blk_idx])
                 blk_idx += 1
                 line_idx += 1
             break
     if not found:
-        raise RuntimeError("The keyword %s is not found in the script." %keyword)
+        return None
     return ret
 
 def get_geometry_in(fname, inlines):
@@ -111,17 +114,21 @@ def get_energy(outlines):
             raise RuntimeError("Final total energy cannot be found in output. Unknown problem.")
     return Etot
 
-def get_force (outlines):
+def get_force (outlines, natoms):
     force = []
-    force_inlines = get_block (outlines, "TOTAL-FORCE (eV/Angstrom)", skip = 4)
+    force_inlines = get_block (outlines, "TOTAL-FORCE (eV/Angstrom)", skip = 4, nlines=np.sum(natoms))
+    if force_inlines is None:
+        raise RuntimeError("TOTAL-FORCE (eV/Angstrom) is not found in running_scf.log. Please check.")
     for line in force_inlines:
         force.append([float(f) for f in line.split()[1:4]])
     force = np.array(force)
     return force
 
 def get_stress(outlines):
     stress = []
-    stress_inlines = get_block(outlines, "TOTAL-STRESS (KBAR)", skip = 3)
+    stress_inlines = get_block(outlines, "TOTAL-STRESS (KBAR)", skip = 3, nlines=3)
+    if stress_inlines is None:
+        return None
     for line in stress_inlines:
         stress.append([float(f) for f in line.split()])
     stress = np.array(stress) * kbar2evperang3
@@ -151,8 +158,10 @@ def get_frame (fname):
     atom_names, natoms, types, coords = get_coords(celldm, cell, geometry_inlines, inlines) 
 
     energy = get_energy(outlines) 
-    force = get_force (outlines) 
-    stress = get_stress(outlines) * np.linalg.det(cell) 
+    force = get_force (outlines, natoms) 
+    stress = get_stress(outlines) 
+    if stress is not None:
+        stress *= np.abs(np.linalg.det(cell)) 
 
     data = {}
     data['atom_names'] = atom_names
@@ -162,7 +171,8 @@ def get_frame (fname):
     data['coords'] = coords[np.newaxis, :, :]
     data['energies'] = np.array(energy)[np.newaxis]
     data['forces'] = force[np.newaxis, :, :]
-    data['virials'] = stress[np.newaxis, :, :]
+    if stress is not None:
+        data['virials'] = stress[np.newaxis, :, :]
     data['orig'] = np.zeros(3)
     # print("atom_names = ", data['atom_names'])
     # print("natoms = ", data['atom_numbs'])

dpdata/plugins/abacus.py

@@ -5,6 +5,7 @@
 
 @Format.register("abacus/scf")
 @Format.register("abacus/pw/scf")
+@Format.register("abacus/lcao/scf")
 class AbacusSCFFormat(Format):
     #@Format.post("rot_lower_triangular")
     def from_labeled_system(self, file_name, **kwargs):

dpdata/qe/scf.py

@@ -70,7 +70,7 @@ def get_coords (lines, cell) :
                 coord.append([float(jj) for jj in ii.split()[1:4]])
                 atom_symbol_list.append(ii.split()[0])
             coord = np.array(coord)
-            coord = np.matmul(coord, cell.T)
+            coord = np.matmul(coord, cell)
     atom_symbol_list = np.array(atom_symbol_list)
     tmp_names, symbol_idx = np.unique(atom_symbol_list, return_index=True)
     atom_types = []

tests/abacus.md/INPUT

@@ -1,30 +1,43 @@
 INPUT_PARAMETERS
-#Parameters	(General)
-suffix          autotest
-pseudo_dir		./
-ntype			1	
-nbands			8
-calculation     md
-read_file_dir   ./
+#Parameters (1.General)
+suffix                  abacus
+calculation             md
+ntype                   2
+nbands                  6
+symmetry                0
 
-#Parameters (Accuracy)
-ecutwfc			20
-niter			20
+#Parameters (2.Iteration)
+ecutwfc                 50
+scf_thr                     1e-2
+scf_nmax                   50
 
-basis_type		pw
-nstep           21 # number of MD/relaxation steps
+#Parameters (3.Basis)
+basis_type              lcao 
+gamma_only              1
 
-stress          1
-stress_thr      1e-6
-force           1
-force_thr_ev    1.0e-3
+#Parameters (4.Smearing)
+smearing_method                gaussian
+smearing_sigma                   0.02
 
-ks_solver       cg
-mixing_type     pulay
-mixing_beta     0.7
+#Parameters (5.Mixing)
+mixing_type             pulay
+mixing_beta             0.4
 
-md_mdtype        1 # 0 for NVE; 1 for NVT with Nose Hoover; 2 for NVT with velocity rescaling
-md_tfirst        10 # temperature, unit: K
-md_dt            1 # unit: fs
-md_rstmd         0 # 1 for restart
-md_dumpmdfred    5 
+#Parameters (6.Deepks)
+#cal_force			1
+#test_force			1
+#deepks_out_labels          1
+#deepks_descriptor_lmax			2
+#deepks_scf	1
+#deepks_model	model.ptg
+
+md_type          1
+md_dt              0.1
+md_tfirst          300
+md_tfreq           0.1
+md_restart           0
+md_dumpfreq          5
+md_nstep               20
+out_level           m
+
+cal_stress 1