merge recent development on devel to master

.github/workflows/test.yml

@@ -0,0 +1,24 @@
+name: Python package
+
+on:
+  - push
+  - pull_request
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: [3.6, 3.7, 3.8]
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@v2
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install dependencies
+      run: pip install . coverage codecov
+    - name: Test
+      run: cd tests && coverage run --source=../dpdata -m unittest && cd .. && coverage combine tests/.coverage && coverage report
+    - run: codecov

README.md

@@ -197,3 +197,13 @@ perturbed_system = dpdata.System('./POSCAR').perturb(pert_num=3,
     atom_pert_style='normal')
 print(perturbed_system.data)
 ```
+
+## replace
+By the following example, Random 8 Hf atoms in the system will be replaced by Zr atoms with the atom postion unchanged.
+```python
+s=dpdata.System('tests/poscars/POSCAR.P42nmc',fmt='vasp/poscar')
+s.replace('Hf', 'Zr', 8)
+s.to_vasp_poscar('POSCAR.P42nmc.replace')
+```
+
+

dpdata/cp2k/output.py

@@ -190,9 +190,12 @@ def handle_single_xyz_frame(self, lines):
                 element_index +=1
                 element_dict[line_list[0]]=[element_index,1]
             atom_types_list.append(element_dict[line_list[0]][0])
-            coords_list.append([float(line_list[1])*AU_TO_ANG,
-                float(line_list[2])*AU_TO_ANG,
-                float(line_list[3])*AU_TO_ANG])
+            # coords_list.append([float(line_list[1])*AU_TO_ANG,
+            #     float(line_list[2])*AU_TO_ANG,
+            #     float(line_list[3])*AU_TO_ANG])
+            coords_list.append([float(line_list[1]),
+                float(line_list[2]),
+                float(line_list[3])])
         atom_names=list(element_dict.keys())
         atom_numbs=[]
         for ii in atom_names:
@@ -210,17 +213,22 @@ def handle_single_xyz_frame(self, lines):
 def get_frames (fname) :
     coord_flag = False
     force_flag = False
+    stress_flag = False
     eV = 2.72113838565563E+01 # hatree to eV
-    angstrom = 5.29177208590000E-01 # Bohrto Angstrom
+    angstrom = 5.29177208590000E-01 # Bohr to Angstrom
+    GPa = 160.21766208 # 1 eV/(Angstrom^3) = 160.21 GPa
     fp = open(fname)
     atom_symbol_list = []
     cell = []
     coord = []
     force = []
+    stress = []
+    cell_count = 0
     coord_count = 0
     for idx, ii in enumerate(fp) :
-        if 'CELL| Vector' in ii :
+        if ('CELL| Vector' in ii) and (cell_count < 3) :
             cell.append(ii.split()[4:7])
+            cell_count += 1
         if 'Atom  Kind  Element' in ii :
             coord_flag = True
             coord_idx = idx
@@ -244,6 +252,18 @@ def get_frames (fname) :
                     force_flag = False
                 else :
                     force.append(ii.split()[3:6])
+        # add reading stress tensor
+        if 'STRESS TENSOR [GPa' in ii :
+            stress_flag = True
+            stress_idx = idx
+        if stress_flag :
+            if (idx > stress_idx + 2):
+                if (ii == '\n') :
+                    stress_flag = False
+                else :
+                    stress.append(ii.split()[1:4])
+
+
     fp.close()
     assert(coord), "cannot find coords"
     assert(energy), "cannot find energies"
@@ -260,10 +280,27 @@ def get_frames (fname) :
     force = np.array(force)
     force = force.astype(np.float)
     force = force[np.newaxis, :, :]
+
+    # virial is not necessary
+    if stress:
+        stress = np.array(stress)
+        stress = stress.astype(np.float)
+        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(cell[0])/GPa
+    else:
+        virial = None
+
+    # force unit conversion, default unit in cp2k is hartree/bohr
     force = force * eV / angstrom
+    # energy unit conversion, default unit in cp2k is hartree
     energy = float(energy) * eV
     energy = np.array(energy)
     energy = energy[np.newaxis]
+
+
+
     tmp_names, symbol_idx = np.unique(atom_symbol_list, return_index=True)
     atom_types = []
     atom_numbs = []
@@ -278,7 +315,7 @@ def get_frames (fname) :
 
     atom_types = np.array(atom_types)
 
-    return list(atom_names), atom_numbs, atom_types, cell, coord, energy, force
+    return list(atom_names), atom_numbs, atom_types, cell, coord, energy, force, virial
 
 
 

dpdata/md/rdf.py

@@ -21,6 +21,15 @@ def rdf(sys,
         Maximal range of rdf calculation
     nbins: int
         Number of bins for rdf calculation
+
+    Returns
+    -------
+    xx: np.array
+        The lattice of r
+    rdf: np.array
+        The value of rdf at r
+    coord: np.array
+        The coordination number up to r
     """
     return compute_rdf(sys['cells'], sys['coords'], sys['atom_types'], 
                        sel_type = sel_type, 
@@ -36,12 +45,16 @@ def compute_rdf(box,
     nframes = box.shape[0]
     xx = None
     all_rdf = []
+    all_cod = []
     for ii in range(nframes):
-        xx, rdf = _compute_rdf_1frame(box[ii], posis[ii], atype, sel_type, max_r, nbins)
+        xx, rdf, cod = _compute_rdf_1frame(box[ii], posis[ii], atype, sel_type, max_r, nbins)
         all_rdf.append(rdf)
+        all_cod.append(cod)
     all_rdf = np.array(all_rdf).reshape([nframes, -1])
+    all_cod = np.array(all_cod).reshape([nframes, -1])
     all_rdf = np.average(all_rdf, axis = 0)
-    return xx, all_rdf
+    all_cod = np.average(all_cod, axis = 0)
+    return xx, all_rdf, all_cod
 
 def _compute_rdf_1frame(box, 
                         posis, 
@@ -65,6 +78,7 @@ def _compute_rdf_1frame(box,
     nlist = ase.neighborlist.NeighborList(max_r, self_interaction=False, bothways=True, primitive=ase.neighborlist.NewPrimitiveNeighborList)
     nlist.update(atoms)
     stat = np.zeros(nbins)
+    stat_acc = np.zeros(nbins)
     hh = max_r / float(nbins)
     for ii in range(natoms) :
         # atom "0"
@@ -90,13 +104,17 @@ def _compute_rdf_1frame(box,
     for ii in sel_type[1]:
         c1 += np.sum(atype == ii)
     rho1 = c1 / np.linalg.det(box)
+    # compute coordination number
+    for ii in range(1, nbins):
+        stat_acc[ii] = stat_acc[ii-1] + stat[ii-1]
+    stat_acc = stat_acc / c0
     # compute rdf
     for ii in range(nbins):
         vol = 4./3. * np.pi * ( ((ii+1)*hh) ** 3 - ((ii)*hh) ** 3 )
         rho = stat[ii] / vol
         stat[ii] = rho / rho1 / c0
     xx = np.arange(0, max_r-1e-12, hh)
-    return xx, stat
+    return xx, stat, stat_acc
 
 if __name__ == '__main__':
     import dpdata

dpdata/md/water.py

@@ -1,5 +1,6 @@
 import numpy as np
 from .pbc import posi_diff
+from .pbc import posi_shift
 
 def compute_bonds (box, 
                    posis,
@@ -179,3 +180,26 @@ def find_ions (atype,
             else :
                 raise RuntimeError("unknow case: numb of O bonded to H > 1")
     return no, noh, noh2, noh3, nh
+
+
+
+def pbc_coords(cells, 
+               coords, 
+               atom_types, 
+               oh_sel = [0, 1],
+               max_roh = 1.3):
+    bonds = compute_bonds(cells, coords, atom_types, oh_sel = oh_sel, max_roh = max_roh, uniq_hbond = True)
+
+    new_coords = np.copy(coords)
+    natoms = len(atom_types)
+    o_type = oh_sel[0]
+    h_type = oh_sel[1]
+    for ii in range(natoms):
+        if atom_types[ii] == o_type:
+            assert(len(bonds[ii]) == 2), 'O has more than 2 bonded Hs, stop'
+            for jj in bonds[ii]:
+                assert(atom_types[jj] == h_type), 'The atom bonded to O is not H, stop'
+                shift = posi_shift(cells, coords[jj], coords[ii])
+                new_coords[jj] = coords[jj] + np.matmul(shift, cells)
+    return new_coords
+

dpdata/qe/scf.py

@@ -64,10 +64,11 @@ def get_coords (lines) :
     return list(atom_names), atom_numbs, atom_types, coord
 
 def get_energy (lines) :
+    energy = None
     for ii in lines :
         if '!    total energy' in ii :
-            return ry2ev * float(ii.split('=')[1].split()[0])
-    return None
+            energy = ry2ev * float(ii.split('=')[1].split()[0])
+    return energy
 
 def get_force (lines) :
     blk = get_block(lines, 'Forces acting on atoms', skip = 1)