merge devel into master

docs/credits.rst

@@ -0,0 +1,4 @@
+Authors
+=======
+
+.. git-shortlog-authors::

docs/index.rst

@@ -13,6 +13,7 @@ Welcome to dpdata's documentation!
    cli
    formats
    api/api
+   credits
 
 .. mdinclude:: ../README.md
 

dpdata/abacus/relax.py

@@ -68,6 +68,7 @@ def get_coords_from_log(loglines,natoms):
             else:
                 cells.append(cells[-1])
 
+            coords[-1] = np.array(coords[-1])
             if direct_coord:
                 coords[-1] = coords[-1].dot(cells[-1])
 

dpdata/abacus/scf.py

@@ -104,15 +104,11 @@ def get_energy(outlines):
             Etot = float(line.split()[1]) # in eV
             break
     if not Etot:
-        not_converge = False
-        for line in outlines:
-            if "convergence has NOT been achieved!" in line:
-                not_converge = True
-                raise RuntimeError("convergence has NOT been achieved in scf!")
-                break
-        if not not_converge:
-            raise RuntimeError("Final total energy cannot be found in output. Unknown problem.")
-    return Etot
+       raise RuntimeError("Final total energy cannot be found in output. Unknown problem.")
+    for line in outlines:
+        if "convergence has NOT been achieved!" in line:
+            return Etot,False
+    return Etot,True
 
 def get_force (outlines, natoms):
     force = []
@@ -156,7 +152,15 @@ def get_frame (fname):
     celldm, cell = get_cell(geometry_inlines) 
     atom_names, natoms, types, coords = get_coords(celldm, cell, geometry_inlines, inlines) 
 
-    energy = get_energy(outlines) 
+    energy,converge = get_energy(outlines) 
+    if not converge:
+        return {'atom_names':atom_names,\
+                'atom_numbs':natoms,\
+                'atom_types':types,\
+                'cells':[],\
+                'coords':[],\
+                'energies':[],\
+                'forces':[]}
     force = get_force (outlines, natoms) 
     stress = get_stress(outlines) 
     if stress is not None:

dpdata/cp2k/output.py

@@ -104,8 +104,8 @@ def handle_single_log_frame(self, lines):
         #  CONSERVED QUANTITY [hartree] =                              -0.279168013085E+04
         energy_pattern_2 = re.compile(r' POTENTIAL ENERGY\[hartree\]\s+=\s+(?P<number>\S+)')
         energy=None
-        cell_length_pattern = re.compile(r' INITIAL CELL LNTHS\[bohr\]\s+=\s+(?P<A>\S+)\s+(?P<B>\S+)\s+(?P<C>\S+)')
-        cell_angle_pattern = re.compile(r' INITIAL CELL ANGLS\[deg\]\s+=\s+(?P<alpha>\S+)\s+(?P<beta>\S+)\s+(?P<gamma>\S+)')
+        cell_length_pattern = re.compile(r' (INITIAL ){0,1}CELL LNTHS\[bohr\]\s+=\s+(?P<A>\S+)\s+(?P<B>\S+)\s+(?P<C>\S+)')
+        cell_angle_pattern = re.compile(r' (INITIAL ){0,1}CELL ANGLS\[deg\]\s+=\s+(?P<alpha>\S+)\s+(?P<beta>\S+)\s+(?P<gamma>\S+)')
         cell_A, cell_B, cell_C = (0,0,0,)
         cell_alpha, cell_beta, cell_gamma=(0,0,0,)
         cell_a_pattern = re.compile(r' CELL\| Vector a \[angstrom\]:\s+(?P<ax>\S+)\s+(?P<ay>\S+)\s+(?P<az>\S+)')

dpdata/driver.py

@@ -147,3 +147,74 @@ def label(self, data: dict) -> dict:
                 labeled_data['energies'] += lb_data ['energies']
                 labeled_data['forces'] += lb_data ['forces']
         return labeled_data
+
+
+class Minimizer(ABC):
+    """The base class for a minimizer plugin. A minimizer can
+    minimize geometry.
+    """
+    __MinimizerPlugin = Plugin()
+
+    @staticmethod
+    def register(key: str) -> Callable:
+        """Register a minimizer plugin. Used as decorators.
+
+        Parameter
+        ---------
+        key: str
+            key of the plugin.
+
+        Returns
+        -------
+        Callable
+            decorator of a class
+
+        Examples
+        --------
+        >>> @Minimizer.register("some_minimizer")
+        ... class SomeMinimizer(Minimizer):
+        ...     pass
+        """
+        return Minimizer.__MinimizerPlugin.register(key)
+
+    @staticmethod
+    def get_minimizer(key: str) -> "Minimizer":
+        """Get a minimizer plugin.
+
+        Parameter
+        ---------
+        key: str
+            key of the plugin.
+
+        Returns
+        -------
+        Minimizer
+            the specific minimizer class
+
+        Raises
+        ------
+        RuntimeError
+            if the requested minimizer is not implemented
+        """
+        try:
+            return Minimizer.__MinimizerPlugin.plugins[key]
+        except KeyError as e:
+            raise RuntimeError('Unknown minimizer: ' + key) from e
+
+    def __init__(self, *args, **kwargs) -> None:
+        """Setup the minimizer."""
+
+    @abstractmethod
+    def minimize(self, data: dict) -> dict:
+        """Minimize the geometry.
+
+        Parameters
+        ----------
+        data : dict
+            data with coordinates and atom types
+
+        Returns
+        -------
+        dict
+            labeled data with minimized coordinates, energies, and forces
+        """

dpdata/fhi_aims/output.py

@@ -1,5 +1,6 @@
 import numpy as np
 import re
+import warnings
 
 latt_patt="\|\s+([0-9]{1,}[.][0-9]*)\s+([0-9]{1,}[.][0-9]*)\s+([0-9]{1,}[.][0-9]*)"
 pos_patt_first="\|\s+[0-9]{1,}[:]\s\w+\s(\w+)(\s.*[-]?[0-9]{1,}[.][0-9]*)(\s+[-]?[0-9]{1,}[.][0-9]*)(\s+[-]?[0-9]{1,}[.][0-9]*)"
@@ -63,7 +64,7 @@ def get_fhi_aims_block(fp) :
             return blk
     return blk
 
-def get_frames (fname, md=True, begin = 0, step = 1) :
+def get_frames (fname, md=True, begin = 0, step = 1, convergence_check=True) :
     fp = open(fname)
     blk = get_fhi_aims_block(fp)
     ret = get_info(blk, type_idx_zero = True)
@@ -78,6 +79,7 @@ def get_frames (fname, md=True, begin = 0, step = 1) :
     all_virials = []    
 
     cc = 0
+    rec_failed = []
     while len(blk) > 0 :
         if debug:
            with open(str(cc),'w') as f:
@@ -87,9 +89,9 @@ def get_frames (fname, md=True, begin = 0, step = 1) :
                 coord, _cell, energy, force, virial, is_converge = analyze_block(blk, first_blk=True, md=md)
             else:
                 coord, _cell, energy, force, virial, is_converge = analyze_block(blk, first_blk=False)
-            if is_converge : 
-                if len(coord) == 0:
-                    break
+            if len(coord) == 0:
+                break
+            if is_converge or not convergence_check: 
                 all_coords.append(coord)
 
                 if _cell:
@@ -101,9 +103,16 @@ def get_frames (fname, md=True, begin = 0, step = 1) :
                 all_forces.append(force)
                 if virial is not None :
                     all_virials.append(virial)
+            if not is_converge:
+                rec_failed.append(cc+1)
+
         blk = get_fhi_aims_block(fp)
         cc += 1
 
+    if len(rec_failed) > 0 :
+        prt = "so they are not collected." if convergence_check else "but they are still collected due to the requirement for ignoring convergence checks."
+        warnings.warn(f"The following structures were unconverged: {rec_failed}; "+prt)
+
     if len(all_virials) == 0 :
         all_virials = None
     else :

dpdata/lammps/dump.py

@@ -5,6 +5,11 @@
 lib_path = os.path.dirname(os.path.realpath(__file__))
 sys.path.append(lib_path)
 import lmp
+import warnings
+class UnwrapWarning(UserWarning):
+    pass
+warnings.simplefilter('once', UnwrapWarning)
+
 
 def _get_block (lines, key) :
     for idx in range(len(lines)) :
@@ -59,16 +64,17 @@ def get_coordtype_and_scalefactor(keys):
     key_su = ['xsu','ysu','zsu'] #scaled and unfolded,sf = lattice parameter
     lmp_coor_type = [key_pc,key_uc,key_s,key_su]
     sf = [0,0,1,1]
+    uw = [0,1,0,1]  # unwraped or not
     for k in range(4):
         if all(i in keys for i in lmp_coor_type[k]):
-            return lmp_coor_type[k],sf[k]
+            return lmp_coor_type[k], sf[k], uw[k]
 
-def safe_get_posi(lines,cell,orig=np.zeros(3)) :
+def safe_get_posi(lines,cell,orig=np.zeros(3), unwrap=False) :
     blk, head = _get_block(lines, 'ATOMS')
     keys = head.split()
     coord_tp_and_sf = get_coordtype_and_scalefactor(keys)
     assert coord_tp_and_sf is not None, 'Dump file does not contain atomic coordinates!'
-    coordtype, sf = coord_tp_and_sf
+    coordtype, sf, uw = coord_tp_and_sf
     id_idx = keys.index('id') - 2
     xidx = keys.index(coordtype[0])-2
     yidx = keys.index(coordtype[1])-2
@@ -81,8 +87,13 @@ def safe_get_posi(lines,cell,orig=np.zeros(3)) :
     posis.sort()
     posis = np.array(posis)[:,1:4]
     if not sf:
-        posis = (posis-orig)@np.linalg.inv(cell)# Convert to scaled coordinates for unscaled coordinates
-    return (posis%1)@cell # Convert scaled coordinates back to Cartesien coordinates 
+        posis = (posis - orig) @ np.linalg.inv(cell)  # Convert to scaled coordinates for unscaled coordinates
+    if uw and unwrap:
+        return posis @ cell  # convert scaled coordinates back to Cartesien coordinates unwrap at the periodic boundaries
+    else:
+        if uw and not unwrap:
+            warnings.warn(message='Your dump file contains unwrapped coordinates, but you did not specify unwrapping (unwrap = True). The default is wrapping at periodic boundaries (unwrap = False).\n',category=UnwrapWarning)
+        return (posis % 1) @ cell  # Convert scaled coordinates back to Cartesien coordinates with wraping at periodic boundary conditions
 
 def get_dumpbox(lines) :
     blk, h = _get_block(lines, 'BOX BOUNDS')
@@ -145,9 +156,9 @@ def load_file(fname, begin = 0, step = 1) :
                 cc += 1
             if cc >= begin and (cc - begin) % step == 0 :
                 buff.append(line)
-
 
-def system_data(lines, type_map = None, type_idx_zero = True) :
+
+def system_data(lines, type_map = None, type_idx_zero = True, unwrap=False) :
     array_lines = split_traj(lines)
     lines = array_lines[0]
     system = {}
@@ -166,12 +177,12 @@ def system_data(lines, type_map = None, type_idx_zero = True) :
     system['cells'] = [np.array(cell)]
     natoms = sum(system['atom_numbs'])
     system['atom_types'] = get_atype(lines, type_idx_zero = type_idx_zero)
-    system['coords'] = [safe_get_posi(lines, cell, np.array(orig))]
+    system['coords'] = [safe_get_posi(lines, cell, np.array(orig), unwrap)]
     for ii in range(1, len(array_lines)) :
         bounds, tilt = get_dumpbox(array_lines[ii])
         orig, cell = dumpbox2box(bounds, tilt)
         system['cells'].append(cell)
-        system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig)))
+        system['coords'].append(safe_get_posi(array_lines[ii], cell, np.array(orig), unwrap))
     system['cells'] = np.array(system['cells'])
     system['coords'] = np.array(system['coords'])
     return system
@@ -181,7 +192,7 @@ def split_traj(dump_lines) :
     marks = []
     for idx,ii in enumerate(dump_lines) :
         if 'ITEM: TIMESTEP' in ii :
-            marks.append(idx) 
+            marks.append(idx)
     if len(marks) == 0 :
         return None
     elif len(marks) == 1 :
@@ -191,9 +202,9 @@ def split_traj(dump_lines) :
         ret = []
         for ii in marks :
             ret.append(dump_lines[ii:ii+block_size])
-        # for ii in range(len(marks)-1): 
+        # for ii in range(len(marks)-1):
         #     assert(marks[ii+1] - marks[ii] == block_size)
-        return ret    
+        return ret
     return None
 
 

dpdata/plugins/3dmol.py

@@ -0,0 +1,44 @@
+from typing import Tuple
+import numpy as np
+
+from dpdata.format import Format
+from dpdata.xyz.xyz import coord_to_xyz
+
+
+@Format.register("3dmol")
+class Py3DMolFormat(Format):
+    """3DMol format.
+
+    To use this format,  py3Dmol should be installed in advance.
+    """
+    def to_system(self,
+                  data: dict,
+                  f_idx: int = 0,
+                  size: Tuple[int] = (300,300),
+                  style: dict = {"stick":{}, "sphere":{"radius":0.4}},
+                  **kwargs):
+        """Show 3D structure of a frame in jupyter.
+
+        Parameters
+        ----------
+        data : dict
+            system data
+        f_idx : int
+            frame index to show
+        size : tuple[int]
+            (width, height) of the widget
+        style : dict
+            style of 3DMol. Read 3DMol documentation for details.
+
+        Examples
+        --------
+        >>> system.to_3dmol()
+        """
+        import py3Dmol
+        types = np.array(data['atom_names'])[data['atom_types']]
+        xyz = coord_to_xyz(data['coords'][f_idx], types)
+        viewer = py3Dmol.view(width=size[0], height=size[1])
+        viewer.addModel(xyz, 'xyz')
+        viewer.setStyle(style.copy())
+        viewer.zoomTo()
+        return viewer

dpdata/plugins/__init__.py

@@ -14,6 +14,9 @@
         importlib.import_module(module_name, PACKAGE_BASE)
 
 # https://setuptools.readthedocs.io/en/latest/userguide/entry_point.html
-eps = metadata.entry_points().get('dpdata.plugins', [])
+try:
+    eps = metadata.entry_points(group='dpdata.plugins')
+except TypeError:
+    eps = metadata.entry_points().get('dpdata.plugins', [])
 for ep in eps:
     plugin = ep.load()

dpdata/plugins/amber.py

@@ -5,7 +5,7 @@
 import dpdata.amber.md
 import dpdata.amber.sqm
 from dpdata.format import Format
-from dpdata.driver import Driver
+from dpdata.driver import Driver, Minimizer
 
 
 @Format.register("amber/md")
@@ -122,3 +122,21 @@ def label(self, data: dict) -> dict:
                             ) from e
                 labeled_system.append(dpdata.LabeledSystem(out_fn, fmt="sqm/out"))
         return labeled_system.data
+
+
+@Minimizer.register("sqm")
+class SQMMinimizer(Minimizer):
+    """SQM minimizer.
+
+    Parameters
+    ----------
+    maxcyc : int, default=1000
+        maximun cycle to minimize
+    """
+    def __init__(self, maxcyc=1000, *args, **kwargs) -> None:
+        assert maxcyc > 0, "maxcyc should be more than 0 to minimize"
+        self.driver = SQMDriver(maxcyc=maxcyc, **kwargs)
+
+    def minimize(self, data: dict) -> dict:
+        # sqm has minimize feature
+        return self.driver.label(data)