merge devel into master

.github/workflows/pub-pypi.yml

@@ -0,0 +1,39 @@
+name: Publish Python distributions to PyPI
+
+on: push
+
+jobs:
+  build-n-publish:
+    if: github.repository_owner == 'deepmodeling'
+    name: Build and publish Python distributions to PyPI
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@master
+
+    - name: Set up Python 3.9
+      uses: actions/setup-python@master
+      with:
+        python-version: 3.9
+
+    - name: Install pypa/build
+      run: >-
+        python -m
+        pip install
+        build
+        --user
+
+    - name: Build a binary wheel and a source tarball
+      run: >-
+        python -m
+        build
+        --sdist
+        --wheel
+        --outdir dist/
+        .
+
+    - name: Publish distribution to PyPI
+      if: startsWith(github.ref, 'refs/tags')
+      uses: pypa/gh-action-pypi-publish@master
+      with:
+        password: ${{ secrets.PYPI_API_TOKEN }}
+

.github/workflows/test.yml

@@ -15,16 +15,15 @@ jobs:
     - uses: actions/checkout@v2
     # set up conda
     - name: Set up Python ${{ matrix.python-version }}
-      uses: conda-incubator/setup-miniconda@v2
+      uses: actions/setup-python@v4
       with:
-        auto-activate-base: true
-        activate-environment: ""
+        python-version: ${{ matrix.python-version }}
     # install rdkit and openbabel
     - name: Install rdkit
-      run: conda create -c conda-forge -n my-rdkit-env python=${{ matrix.python-version }} rdkit openbabel;
+      run: python -m pip install rdkit openbabel-wheel
     - name: Install dependencies
-      run: source $CONDA/bin/activate my-rdkit-env && pip install .[amber,ase,pymatgen] coverage codecov
+      run: python -m pip install .[amber,ase,pymatgen] coverage codecov
     - name: Test
-      run: source $CONDA/bin/activate my-rdkit-env && cd tests && coverage run --source=../dpdata -m unittest && cd .. && coverage combine tests/.coverage && coverage report
+      run: cd tests && coverage run --source=../dpdata -m unittest && cd .. && coverage combine tests/.coverage && coverage report
     - name: Run codecov
-      run: source $CONDA/bin/activate my-rdkit-env && codecov
+      run: codecov

README.md

@@ -82,6 +82,7 @@ The `System` or `LabeledSystem` can be constructed from the following file forma
 | Gromacs | gro         | True         | False   | System        | 'gromacs/gro'        |
 | ABACUS  | STRU        | False        | True    | LabeledSystem | 'abacus/scf'         |
 | ABACUS  | cif         | True         | True    | LabeledSystem | 'abacus/md'          |
+| ABACUS  | STRU        | True         | True    | LabeledSystem | 'abacus/relax'       |
 | ase     | structure   | True         | True    | MultiSystems  | 'ase/structure'      |
 
 

docs/cli.rst

@@ -0,0 +1,7 @@
+Command line interface
+======================
+
+.. argparse::
+   :module: dpdata.cli
+   :func: dpdata_parser
+   :prog: dpdata

docs/conf.py

@@ -48,6 +48,7 @@
     'sphinx.ext.intersphinx',
     'numpydoc',
     'm2r2',
+    'sphinxarg.ext',
 ]
 
 # Add any paths that contain templates here, relative to this directory.

docs/index.rst

@@ -10,6 +10,7 @@ Welcome to dpdata's documentation!
    :maxdepth: 2
    :caption: Contents:
 
+   cli
    formats
    api/api
 

dpdata/abacus/relax.py

@@ -0,0 +1,139 @@
+import os,sys
+import numpy as np
+from .scf import bohr2ang, kbar2evperang3, get_geometry_in, get_cell, get_coords
+
+# Read in geometries from an ABACUS RELAX(CELL-RELAX) trajectory in OUT.XXXX/runnning_relax/cell-relax.log. 
+
+def get_log_file(fname, inlines):
+    suffix = "ABACUS"
+    calculation = "scf"
+    for line in inlines:
+        if "suffix" in line and "suffix"==line.split()[0]:
+           suffix = line.split()[1]
+        elif "calculation" in line and "calculation" == line.split()[0]:
+            calculation = line.split()[1]
+    logf = os.path.join(fname, "OUT.%s/running_%s.log"%(suffix,calculation))
+    return logf
+
+def get_coords_from_log(loglines,natoms):
+    '''
+    NOTICE: unit of coords and cells is Angstrom
+    '''
+    natoms_log = 0
+    for line in loglines:
+        if line[13:41] == "number of atom for this type":
+            natoms_log += int(line.split()[-1])
+
+    assert(natoms_log>0 and natoms_log == natoms),"ERROR: detected atom number in log file is %d" % natoms
+
+    energy = []
+    cells = []
+    coords = []
+    force = []
+    stress = []
+
+    for i in range(len(loglines)):
+        line = loglines[i]
+        if line[18:41] == "lattice constant (Bohr)": 
+            a0 = float(line.split()[-1])
+        elif len(loglines[i].split()) >=2 and loglines[i].split()[1] == 'COORDINATES':
+            coords.append([])
+            direct_coord = False
+            if loglines[i].split()[0] == 'DIRECT':
+                direct_coord = True
+                for k in range(2,2+natoms):
+                    coords[-1].append(list(map(lambda x: float(x),loglines[i+k].split()[1:4])))
+            elif loglines[i].split()[0] == 'CARTESIAN':
+                for k in range(2,2+natoms):
+                    coords[-1].append(list(map(lambda x: float(x)*a0,loglines[i+k].split()[1:4])))  
+            else:
+                assert(False),"Unrecongnized coordinate type, %s, line:%d" % (loglines[i].split()[0],i)
+
+            converg = True
+            for j in range(i):
+                if loglines[i-j-1][1:36] == 'Ion relaxation is not converged yet':
+                    converg = False
+                    break
+                elif loglines[i-j-1][1:29] == 'Ion relaxation is converged!':
+                    converg = True
+                    break
+
+            if converg:
+                for j in range(i+1,len(loglines)):
+                    if loglines[j][1:56] == "Lattice vectors: (Cartesian coordinate: in unit of a_0)":
+                        cells.append([])
+                        for k in range(1,4):
+                            cells[-1].append(list(map(lambda x:float(x)*a0,loglines[j+k].split()[0:3]))) 
+                        break
+            else:
+                cells.append(cells[-1])
+
+            if direct_coord:
+                coords[-1] = coords[-1].dot(cells[-1])
+
+        elif line[4:15] == "TOTAL-FORCE":
+            force.append([])
+            for j in range(5,5+natoms):
+                force[-1].append(list(map(lambda x:float(x),loglines[i+j].split()[1:4])))
+        elif line[1:13] == "TOTAL-STRESS":
+            stress.append([])
+            for j in range(4,7):
+                stress[-1].append(list(map(lambda x:float(x),loglines[i+j].split()[0:3])))
+        elif line[1:14] == "final etot is":
+            energy.append(float(line.split()[-2]))
+
+    assert(len(cells) == len(coords) or len(cells)+1 == len(coords)),"ERROR: detected %d coordinates and %d cells" % (len(coords),len(cells))
+    if len(cells)+1 == len(coords): del(coords[-1])
+
+    energy = np.array(energy)
+    cells = np.array(cells)
+    coords = np.array(coords)
+    stress = np.array(stress)
+    force = np.array(force)
+
+    cells *= bohr2ang
+    coords *= bohr2ang
+
+    virial = np.zeros([len(cells), 3, 3])
+    for i in range(len(cells)):
+        volume = np.linalg.det(cells[i, :, :].reshape([3, 3]))
+        virial[i] = stress[i] * kbar2evperang3 * volume
+
+    return energy,cells,coords,force,stress,virial
+
+def get_frame (fname):
+    if type(fname) == str:
+        # if the input parameter is only one string, it is assumed that it is the 
+        # base directory containing INPUT file;
+        path_in = os.path.join(fname, "INPUT")
+    else:
+        raise RuntimeError('invalid input')    
+    with open(path_in, 'r') as fp:
+        inlines = fp.read().split('\n')
+    geometry_path_in = get_geometry_in(fname, inlines) # base dir of STRU
+    with open(geometry_path_in, 'r') as fp:
+        geometry_inlines = fp.read().split('\n')
+    celldm, cell = get_cell(geometry_inlines) 
+    atom_names, natoms, types, coord_tmp = get_coords(celldm, cell, geometry_inlines, inlines) 
+
+    logf = get_log_file(fname, inlines) 
+    assert(os.path.isfile(logf)),"Error: can not find %s" % logf
+    with open(logf) as f1: lines = f1.readlines()
+
+    atomnumber = 0
+    for i in natoms: atomnumber += i
+    energy,cells,coords,force,stress,virial = get_coords_from_log(lines,atomnumber)
+
+    data = {}
+    data['atom_names'] = atom_names
+    data['atom_numbs'] = natoms
+    data['atom_types'] = types
+    data['cells'] = cells
+    data['coords'] = coords
+    data['energies'] = energy
+    data['forces'] = force
+    data['virials'] = virial
+    data['stress'] = stress
+    data['orig'] = np.zeros(3)
+
+    return data

dpdata/abacus/scf.py

@@ -33,7 +33,7 @@ def get_block (lines, keyword, skip = 0, nlines = None):
 def get_geometry_in(fname, inlines):
     geometry_path_in = os.path.join(fname, "STRU")
     for line in inlines:
-        if "atom_file" in line and "atom_file"==line.split()[0]:
+        if "stru_file" in line and "stru_file"==line.split()[0]:
            atom_file = line.split()[1]
            geometry_path_in = os.path.join(fname, atom_file)
            break
@@ -148,7 +148,6 @@ def get_frame (fname):
 
     geometry_path_in = get_geometry_in(fname, inlines) 
     path_out = get_path_out(fname, inlines) 
-
     with open(geometry_path_in, 'r') as fp:
         geometry_inlines = fp.read().split('\n')
     with open(path_out, 'r') as fp:
@@ -233,7 +232,7 @@ def make_unlabeled_stru(data, frame_idx, pp_file=None, numerical_orbital=None, n
     for iele in range(len(data['atom_names'])):
         out += data['atom_names'][iele] + " "
         if mass is not None:
-            out += "%d "%mass[iele]
+            out += "%.3f "%mass[iele]
         else:
             out += "1 "
         if pp_file is not None:
@@ -280,4 +279,4 @@ def make_unlabeled_stru(data, frame_idx, pp_file=None, numerical_orbital=None, n
 
 #if __name__ == "__main__":
 #    path = "/home/lrx/work/12_ABACUS_dpgen_interface/dpdata/dpdata/tests/abacus.scf"
-#    data = get_frame(path)
+#    data = get_frame(path)

dpdata/amber/sqm.py

@@ -94,12 +94,12 @@ def make_sqm_in(data, fname=None, frame_idx=0, **kwargs):
     ret += "     verbosity=4\n"
     ret += " /\n"
     for ii in range(len(data['atom_types'])):
-        ret += "{:>4s}{:>6s}{:>14s}{:>14s}{:>14s}\n".format(
+        ret += "{:>4s}{:>6s}{:>16s}{:>16s}{:>16s}\n".format(
             str(atomic_numbers[ii]),
             str(symbols[ii]),
-            f"{data['coords'][frame_idx][ii, 0]:.4f}",
-            f"{data['coords'][frame_idx][ii, 1]:.4f}",
-            f"{data['coords'][frame_idx][ii, 2]:.4f}"
+            f"{data['coords'][frame_idx][ii, 0]:.6f}",
+            f"{data['coords'][frame_idx][ii, 1]:.6f}",
+            f"{data['coords'][frame_idx][ii, 2]:.6f}"
         )
     if fname is not None:
         with open(fname, 'w') as fp:

dpdata/cli.py

@@ -5,14 +5,13 @@
 from .system import System, LabeledSystem, MultiSystems
 
 
-def dpdata_cli():
-    """dpdata cli.
+def dpdata_parser() -> argparse.ArgumentParser:
+    """Returns dpdata cli parser.
 
-    Examples
-    --------
-    .. code-block:: bash
-
-        $ dpdata -iposcar POSCAR -odeepmd/npy -O data -n
+    Returns
+    -------
+    argparse.ArgumentParser
+        dpdata cli parser
     """
     parser = argparse.ArgumentParser(
         description="dpdata: Manipulating multiple atomic simulation data formats",
@@ -28,7 +27,19 @@ def dpdata_cli():
     parser.add_argument("--type-map", "-t", type=str, nargs="+", help="type map")
 
     parser.add_argument('--version', action='version', version='dpdata v%s' % __version__)
+    return parser
 
+
+def dpdata_cli():
+    """dpdata cli.
+
+    Examples
+    --------
+    .. code-block:: bash
+
+        $ dpdata -iposcar POSCAR -odeepmd/npy -O data -n
+    """
+    parser = dpdata_parser()
     parsed_args = parser.parse_args()
     convert(**vars(parsed_args))
 

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):