Merge devel into master (#385)

Author: amcadmus

.github/workflows/test.yml

@@ -22,8 +22,13 @@ jobs:
     - name: Install rdkit
       run: python -m pip install rdkit openbabel-wheel
     - name: Install dependencies
-      run: python -m pip install .[amber,ase,pymatgen] coverage codecov
+      run: python -m pip install .[amber,ase,pymatgen] coverage
     - name: Test
       run: cd tests && coverage run --source=../dpdata -m unittest && cd .. && coverage combine tests/.coverage && coverage report
     - name: Run codecov
-      run: codecov
+      uses: codecov/codecov-action@v3
+  pass:
+    needs: [build]
+    runs-on: ubuntu-latest
+    steps:
+      - run: echo "All jobs passed"

README.md

@@ -259,11 +259,8 @@ If a valence of 3 is detected on carbon, the formal charge will be assigned to -
 
 # Plugins
 
-One can follow [a simple example](plugin_example/) to add their own format by creating and installing plugins. It's critical to add the [Format](dpdata/format.py) class to `entry_points['dpdata.plugins']` in `setup.py`:
-```py
-    entry_points={
-        'dpdata.plugins': [
-            'random=dpdata_random:RandomFormat'
-        ]
-    },
+One can follow [a simple example](plugin_example/) to add their own format by creating and installing plugins. It's critical to add the [Format](dpdata/format.py) class to `entry_points['dpdata.plugins']` in [`pyproject.toml`](plugin_example/pyproject.toml):
+```toml
+[project.entry-points.'dpdata.plugins']
+random = "dpdata_random:RandomFormat"
 ```

dpdata/abacus/md.py

@@ -3,6 +3,7 @@
 import numpy as np
 from .scf import ry2ev, bohr2ang, kbar2evperang3, get_block, get_geometry_in, get_cell, get_coords
 import re
+import warnings
 
 # Read in geometries from an ABACUS MD trajectory.
 # The atomic coordinates are read in from generated files in OUT.XXXX.
@@ -80,6 +81,10 @@ def get_energy(outlines, ndump, dump_freq):
             if nenergy%dump_freq == 0:
                 energy.append(float(line.split()[-2]))
             nenergy+=1
+        elif "!! convergence has not been achieved" in line:
+            if nenergy%dump_freq == 0:
+                energy.append(np.nan)
+            nenergy+=1
     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
@@ -113,6 +118,20 @@ def get_frame (fname):
     with open(os.path.join(path_out, "running_md.log"), 'r') as fp:
         outlines = fp.read().split('\n')
     energy = get_energy(outlines, ndump, dump_freq)
+
+    unconv_stru = ''
+    for i,iene in enumerate(energy):
+        if np.isnan(iene):
+            coords = np.delete(coords,i-ndump,axis=0)
+            cells = np.delete(cells,i-ndump,axis=0)
+            force = np.delete(force,i-ndump,axis=0)
+            stress = np.delete(stress,i-ndump,axis=0)
+            energy = np.delete(energy,i-ndump,axis=0) 
+            unconv_stru += "%d " % i
+    ndump = len(energy)
+    if unconv_stru != '':
+        warnings.warn(f"Structure %s are unconverged and not collected!" % unconv_stru)  
+
     for iframe in range(ndump):
         stress[iframe] *= np.linalg.det(cells[iframe, :, :].reshape([3, 3]))
     if np.sum(np.abs(stress[0])) < 1e-10:

dpdata/cli.py

@@ -1,5 +1,6 @@
 """Command line interface for dpdata."""
 import argparse
+from typing import Optional
 
 from . import __version__
 from .system import System, LabeledSystem, MultiSystems
@@ -47,11 +48,11 @@ def dpdata_cli():
 def convert(*,
             from_file: str,
             from_format: str = "auto",
-            to_file: str = None,
-            to_format: str = None,
+            to_file: Optional[str] = None,
+            to_format: Optional[str] = None,
             no_labeled: bool = False,
             multi: bool = False,
-            type_map: list = None,
+            type_map: Optional[list] = None,
             **kwargs):
     """Convert files from one format to another one.
     

dpdata/deepmd/hdf5.py

@@ -1,5 +1,5 @@
 """Utils for deepmd/hdf5 format."""
-from typing import Union
+from typing import Optional, Union
 
 import h5py
 import numpy as np
@@ -11,7 +11,7 @@
 
 def to_system_data(f: Union[h5py.File, h5py.Group],
                    folder: str, 
-                   type_map: list = None,
+                   type_map: Optional[list] = None,
                    labels: bool = True) :
     """Load a HDF5 file.
 

dpdata/gaussian/gjf.py

@@ -3,7 +3,7 @@
 # under LGPL 3.0 license
 """Generate Gaussian input file."""
 
-from typing import List, Tuple, Union
+from typing import Optional, List, Tuple, Union
 import uuid
 import itertools
 import warnings
@@ -107,8 +107,8 @@ def make_gaussian_input(
         multiplicity: Union[str ,int] = "auto",
         charge: int = 0,
         fragment_guesses: bool = False,
-        basis_set: str = None,
-        keywords_high_multiplicity: str = None,
+        basis_set: Optional[str] = None,
+        keywords_high_multiplicity: Optional[str] = None,
         nproc: int = 1,
         ) -> str:
     """Make gaussian input file.

dpdata/lammps/dump.py

@@ -31,9 +31,11 @@ def get_atype(lines, type_idx_zero = False) :
     tidx = keys.index('type') - 2
     atype = []
     for ii in blk :
-        atype.append([int(ii.split()[id_idx]), int(ii.split()[tidx])])
+        atype.append([int(ii.split()[tidx]), int(ii.split()[id_idx])])
+    # sort with type id
     atype.sort()
     atype = np.array(atype, dtype = int)    
+    atype = atype[:, ::-1]
     if type_idx_zero :
         return atype[:,1] - 1
     else :
@@ -76,16 +78,17 @@ def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) :
     assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!'
     coordtype, sf, uw = coord_tp_and_sf
     id_idx = keys.index('id') - 2
+    tidx = keys.index('type') - 2
     xidx = keys.index(coordtype[0])-2
     yidx = keys.index(coordtype[1])-2
     zidx = keys.index(coordtype[2])-2
     sel = (xidx, yidx, zidx)
     posis = []
     for ii in blk :
         words = ii.split()
-        posis.append([float(words[id_idx]), float(words[xidx]), float(words[yidx]), float(words[zidx])])
+        posis.append([float(words[tidx]), float(words[id_idx]), float(words[xidx]), float(words[yidx]), float(words[zidx])])
     posis.sort()
-    posis = np.array(posis)[:,1:4]
+    posis = np.array(posis)[:,2:5]
     if not sf:
         posis = (posis - orig) @ np.linalg.inv(cell)  # Convert to scaled coordinates for unscaled coordinates
     if uw and unwrap:

dpdata/plugins/ase.py

@@ -1,4 +1,4 @@
-from typing import TYPE_CHECKING, Type
+from typing import Optional, TYPE_CHECKING, Type
 from dpdata.driver import Driver, Minimizer
 from dpdata.format import Format
 import numpy as np
@@ -93,7 +93,7 @@ def from_labeled_system(self, atoms: "ase.Atoms", **kwargs) -> dict:
             info_dict['virials'] = virials
         return info_dict
 
-    def from_multi_systems(self, file_name: str, begin: int=None, end: int=None, step: int=None, ase_fmt: str=None, **kwargs) -> "ase.Atoms":
+    def from_multi_systems(self, file_name: str, begin: Optional[int] = None, end: Optional[int] = None, step: Optional[int] = None, ase_fmt: Optional[str] = None, **kwargs) -> "ase.Atoms":
         """Convert a ASE supported file to ASE Atoms.
 
         It will finally be converted to MultiSystems.
@@ -221,13 +221,16 @@ class ASEMinimizer(Minimizer):
         ase optimizer class
     fmax : float, optional, default=5e-3
         force convergence criterion
+    max_steps : int, optional
+        max steps to optimize
     optimizer_kwargs : dict, optional
         other parameters for optimizer
     """
     def __init__(self,
                  driver: Driver,
-                 optimizer: Type["Optimizer"] = None,
+                 optimizer: Optional[Type["Optimizer"]] = None,
                  fmax: float = 5e-3,
+                 max_steps: Optional[int] = None,
                  optimizer_kwargs: dict = {}) -> None:
         self.calculator = driver.ase_calculator
         if optimizer is None:
@@ -240,6 +243,7 @@ def __init__(self,
             **optimizer_kwargs.copy(),
         }
         self.fmax = fmax
+        self.max_steps = max_steps
 
     def minimize(self, data: dict) -> dict:
         """Minimize the geometry.
@@ -261,7 +265,7 @@ def minimize(self, data: dict) -> dict:
         for atoms in structures:
             atoms.calc = self.calculator
             dyn = self.optimizer(atoms, **self.optimizer_kwargs)
-            dyn.run(fmax=self.fmax)
+            dyn.run(fmax=self.fmax, steps=self.max_steps)
             ls = dpdata.LabeledSystem(atoms, fmt="ase/structure", type_map=data['atom_names'])
             labeled_system.append(ls)
-        return labeled_system.data
+        return labeled_system.data

dpdata/plugins/deepmd.py

@@ -1,4 +1,4 @@
-from typing import Union, List
+from typing import Optional, Union, List
 
 import dpdata
 import dpdata.deepmd.raw
@@ -108,7 +108,7 @@ def _from_system(self, file_name: Union[str, h5py.Group, h5py.File], type_map: L
 
     def from_system(self,
                     file_name: Union[str, h5py.Group, h5py.File],
-                    type_map: List[str]=None,
+                    type_map: Optional[List[str]] = None,
                     **kwargs) -> dict:
         """Convert HDF5 file to System data.
 
@@ -134,7 +134,7 @@ def from_system(self,
 
     def from_labeled_system(self,
                             file_name: Union[str, h5py.Group, h5py.File],
-                            type_map: List[str]=None,
+                            type_map: Optional[List[str]] = None,
                             **kwargs) -> dict:
         """Convert HDF5 file to LabeledSystem data.
 

dpdata/system.py

@@ -810,6 +810,11 @@ def perturb(self,
         perturbed_system : System
             The perturbed_system. It contains `pert_num` * frame_num of the input system frames.
         """
+        if type(self) is not dpdata.System:
+            raise RuntimeError(
+                f'Using method perturb() of an instance of {type(self)}. '
+                f'Must use method perturb() of the instance of class dpdata.System.'
+            )
         perturbed_system = System()
         nframes = self.get_nframes()
         for ii in range(nframes):
@@ -1125,11 +1130,6 @@ def affine_map_fv(self, trans, f_idx) :
         if self.has_virial():
             self.data['virials'][f_idx] = np.matmul(trans.T, np.matmul(self.data['virials'][f_idx], trans))
 
-    @post_funcs.register("rot_lower_triangular")
-    def rot_lower_triangular(self) :
-        for ii in range(self.get_nframes()) :
-            self.rot_frame_lower_triangular(ii)
-
     def rot_frame_lower_triangular(self, f_idx = 0) :
         trans = System.rot_frame_lower_triangular(self, f_idx = f_idx)
         self.affine_map_fv(trans, f_idx = f_idx)