Merge branch 'develop' of https://github.com/abinit/abipy into memory…

…_units

abipy/abilab.py

@@ -423,7 +423,7 @@ def cmp_version(this: str, other: str, op: str = ">=") -> bool:
         Compare two version strings with the given operator ``op``
         >>> assert cmp_version("1.1.1", "1.1.0") and not cmp_version("1.1.1", "1.1.0", op="==")
         """
-        from pkg_resources import parse_version
+        from packaging.version import parse as parse_version
         from monty.operator import operator_from_str
         op = operator_from_str(op)
         return op(parse_version(this), parse_version(other))

abipy/abio/abivar_database/variables_optic.py

@@ -237,7 +237,7 @@
 This parameter provides a fixed shift to all the conduction bands. As
 LDA/GGA are known to underestimate the band-gap by a significant amount in
 some cases, in order to obtain a reasonable optical spectrum and make a realistic
-comparison with experiments one needs to correct for this.
+comparison with experiments one needs to correct for this [[cite:Levine1989]].
 The scissors shift is normally chosen to be the difference between the experimental and
 theoretical band-gap, and simply shifts the conduction bands. Alternatively, one may
 determine the self energy using the [[tutorial:gw1|GW approach]], in which case

abipy/abio/inputs.py

@@ -472,7 +472,7 @@ def enforce_znucl_and_typat(self, znucl, typat):
         if len(typat) != len(self.structure):
             raise ValueError("typat contains %d entries while it should be natom: %d" % (len(typat), len(self.structure)))
 
-        ntypat = self.structure.ntypesp
+        ntypat = self.structure.n_elems
         if len(znucl) != ntypat:
             raise ValueError("znucl contains %d entries while it should be ntypat: %d" % (len(znucl), ntypat))
 
@@ -814,7 +814,7 @@ def escape(text):
                 for name in names:
                     value = vars[name]
                     if mnemonics and value is not None:
-                        print(f"{name=}")
+                        #print(f"{name=}")
                         app(escape("#### <" + var_database[name].mnemonics + ">"))
 
                     # Build variable, convert to string and append it
@@ -4053,6 +4053,7 @@ def to_string(self, sortmode=None, mode="text", verbose=0, files_file=False) ->
             #if mode == "html": vname = root + "#%s" % vname
             if mode == "html": vname = var_database[vname].html_link(label=vname)
             value = format_string_abivars(vname, value, "anaddb")
+            #print("vname:", vname, "value:", value)
             app(str(InputVariable(vname, value)))
 
         return "\n".join(lines) if mode == "text" else "\n".join(lines).replace("\n", "<br>")

abipy/core/kpoints.py

@@ -101,7 +101,6 @@ def issamek(k1, k2, atol=None):
     """
     k1 = np.asarray(k1)
     k2 = np.asarray(k2)
-    #if k1.shape != k2.shape:
 
     return is_integer(k1 - k2, atol=atol)
 
@@ -409,32 +408,35 @@ def map_kpoints(other_kpoints, other_lattice, ref_lattice, ref_kpoints, ref_symr
 def kpoints_indices(frac_coords, ngkpt, check_mesh=0) -> np.ndarray:
     """
     This function is used when we need to insert k-dependent quantities in a (nx, ny, nz) array.
-    It compute the indices of the k-points assuming these points belong to
-    a mesh with ngkpt divisions.
+    It computes the indices of the k-points assuming these points belong to a mesh with ngkpt divisions.
 
     Args:
         frac_coords
         ngkpt:
-        check_mesh
+        check_mesh:
     """
+
     # Transforms kpt in its corresponding reduced number in the interval [0,1[
     k_indices = [np.round((kpt % 1) * ngkpt) for kpt in frac_coords]
-
     k_indices = np.array(k_indices, dtype=int)
-    #k_indices = np.array(list(map(tuple, k_indices)))
-    #k_indices = list(map(tuple, k_indices))
 
+    # Debug secction.
     if check_mesh:
         print(f"kpoints_indices: Testing whether k-points belong to the {ngkpt =} mesh")
         ierr = 0
         for kpt, inds in zip(frac_coords, k_indices):
-            #print("kpt:", kpt, "inds:", inds)
+            if check_mesh > 1: print("kpt:", kpt, "inds:", inds)
             same_k = np.array((inds[0]/ngkpt[0], inds[1]/ngkpt[1], inds[2]/ngkpt[2]))
             if not issamek(kpt, same_k):
                 ierr += 1; print(kpt, "-->", same_k)
         if ierr:
             raise ValueError("Wrong mapping")
-        print("Check succesful!")
+
+        #for kpt, inds in zip(frac_coords, k_indices):
+        #    if np.any(inds >= ngkpt):
+        #        raise ValueError(f"inds >= nkgpt for {kpt=}, {np.round(kpt % 1)=} {inds=})")
+
+        print("Check succesfull!")
 
     return k_indices
 
@@ -526,6 +528,7 @@ def kpath_from_bounds_and_ndivsm(bounds, ndivsm, structure):
         for j in range(ndivs[i]):
             p = bounds[i] + j * (bounds[i + 1] - bounds[i]) / ndivs[i]
             path.append(p)
+
     path.append(bounds[-1])
 
     return np.array(path)
@@ -1400,6 +1403,8 @@ def lines(self):
             for line in self.lines:
                 vals_on_line = eigens[spin, line, band]
         """
+        if len(self) < 2:
+            return tuple()
         prev = self.versors[0]
         lines = [[0]]
 

abipy/core/skw.py

@@ -890,8 +890,8 @@ def __init__(self, lpratio, kpts, eigens, fermie, nelect, cell, symrel, has_timr
         cprint("FIT vs input data: Mean Absolute Error= %.3e (meV)" % mae, color="red" if warn else "green")
         if warn:
             # Issue warning if error too large.
-            cprint("Large error in SKW interpolation!", "red")
-            cprint("MAE:", mae, "[meV]", "red")
+            cprint("Large error in SKW interpolation!", color="red")
+            cprint(f"MAE: {mae} [meV]", color="red")
 
         self.mae = mae
 

abipy/core/structure.py

@@ -1272,6 +1272,50 @@ def print_neighbors(self, radius=2.0) -> None:
                 print("\t", repr(s), " at distance", dist)
             print("")
 
+    @lazy_property
+    def has_zero_dynamical_quadrupoles(self):
+        """
+        Dynamical quadrupoles are nonzero in all noncentrosymmetric crystals,
+        but also in centrosymmetric ones if one or more atoms are placed at noncentrosymmetric sites.
+        """
+        def create_image(s1, s2):
+            """
+            Creates the image of s2 through s1
+            This image is a fictitious atom in the structure
+            """
+            image = PeriodicSite.from_dict(s2.as_dict())
+            image.coords = s1.coords - (s2.coords - s1.coords)
+            image.frac_coords = s1.frac_coords - (s2.frac_coords - s1.frac_coords)
+
+            return image
+
+        def check_image(structure, site):
+            """
+            Checks if a fictitious site is an image of a site of the structure
+            """
+            for site in structure.sites:
+                if site.is_periodic_image(site):
+                    return True
+
+            return False
+
+        # If the centrosymmetry is broken at a given atomic site of the given structure,
+        # returns False. Else, return True
+        sites = self.sites
+
+        for s1 in sites:
+            for s2 in sites:
+                # Do not take s1 = s2 into account
+                if s2 != s1:
+                    # Create the image of s2 through s1 (fictitious atom)
+                    image = create_image(s1, s2)
+                    is_image = check_image(self, image)
+                    if not is_image:
+                        return False
+
+        return True
+
+
     @lazy_property
     def hsym_kpath(self):
         """
@@ -1892,8 +1936,7 @@ def get_smallest_supercell(self, qpoint, max_supercell):
         Returns: the scaling matrix of the supercell
         """
         if np.allclose(qpoint, 0.0):
-            scale_matrix = np.eye(3, 3)
-            return scale_matrix
+            return np.eye(3, 3)
 
         l = max_supercell
 
@@ -1966,11 +2009,12 @@ def get_smallest_supercell(self, qpoint, max_supercell):
             raise ValueError('max_supercell is not large enough for this q-point')
 
         # Fortran 2 python!!!
-        return scale_matrix.T
+        return scale_matrix.T.copy()
 
-    def make_doped_supercells(self,scaling_matrix,replaced_atom, dopant_atom):
+    def make_doped_supercells(self, scaling_matrix, replaced_atom, dopant_atom):
         """
         Returns a list doped supercell structures, one for each non-equivalent site of the replaced atom.
+
         Args:
             scaling_matrix: A scaling matrix for transforming the lattice vectors.
                 Has to be all integers. Several options are possible:
@@ -1985,18 +2029,18 @@ def make_doped_supercells(self,scaling_matrix,replaced_atom, dopant_atom):
             dopant_atom : Symbol of the dopant_atom (ex: 'Eu')
         """
         ### list of positions of non-equivalent sites for the replaced atom. ###
-        irred=self.spget_equivalent_atoms().eqmap # mapping from inequivalent sites to atoms sites
-        positions=self.get_symbol2indices()[replaced_atom] # get indices of the replaced atom
+        irred = self.spget_equivalent_atoms().eqmap # mapping from inequivalent sites to atoms sites
+        positions = self.get_symbol2indices()[replaced_atom] # get indices of the replaced atom
 
-        index_non_eq_sites=[]
+        index_non_eq_sites = []
         for pos in positions:
             if len(irred[pos]) != 0:
                  index_non_eq_sites.append(irred[pos][0])
 
-        doped_supercell=self.copy()
+        doped_supercell = self.copy()
         doped_supercell.make_supercell(scaling_matrix)
 
-        doped_structure_list=[]
+        doped_structure_list = []
 
         for index in index_non_eq_sites:
             final_structure=doped_supercell.copy()
@@ -2005,8 +2049,6 @@ def make_doped_supercells(self,scaling_matrix,replaced_atom, dopant_atom):
 
         return doped_structure_list
 
-
-
     def get_trans_vect(self, scale_matrix):
         """
         Returns the translation vectors for a given scale matrix.
@@ -2038,8 +2080,7 @@ def range_vec(i):
 
         # find the translation vectors (in terms of the initial lattice vectors)
         # that are inside the unit cell defined by the scale matrix
-        # we're using a slightly offset interval from 0 to 1 to avoid numerical
-        # precision issues
+        # we're using a slightly offset interval from 0 to 1 to avoid numerical precision issues
         #print(scale_matrix)
         inv_matrix = np.linalg.inv(scale_matrix)
 
@@ -2106,8 +2147,8 @@ def write_vib_file(self, xyz_file, qpoint, displ, do_real=True, frac_coords=True
     def frozen_2phonon(self, qpoint, displ1, displ2, eta=1, frac_coords=False, scale_matrix=None, max_supercell=None):
         """
         Creates the supercell needed for a given qpoint and adds the displacements.
-        The displacements are normalized so that the largest atomic displacement will correspond to the
-        value of eta in Angstrom.
+        The displacements are normalized so that the largest atomic displacement will correspond
+        to the value of eta in Angstrom.
 
         Args:
             qpoint: q vector in reduced coordinate in reciprocal space.
@@ -2182,15 +2223,14 @@ def frozen_phonon(self, qpoint, displ, eta=1, frac_coords=False, scale_matrix=No
         value of eta in Angstrom.
 
         Args:
-            qpoint: q vector in reduced coordinate in reciprocal space.
+            qpoint: q-vector in reduced coordinate in reciprocal space.
             displ: displacement in real space of the atoms.
-            eta: pre-factor multiplying the displacement. Gives the value in Angstrom of the
-                largest displacement.
+            eta: pre-factor multiplying the displacement. Gives the value in Angstrom of the largest displacement.
             frac_coords: whether the displacements are given in fractional or cartesian coordinates
             scale_matrix: the scaling matrix of the supercell. If None a scaling matrix suitable for
                 the qpoint will be determined.
             max_supercell: mandatory if scale_matrix is None, ignored otherwise. Defines the largest
-                supercell in the search for a scaling matrix suitable for the q point.
+                supercell in the search for a scaling matrix suitable for the input q-point.
 
         Returns:
             A namedtuple with a Structure with the displaced atoms, a numpy array containing the
@@ -2199,35 +2239,34 @@ def frozen_phonon(self, qpoint, displ, eta=1, frac_coords=False, scale_matrix=No
 
         if scale_matrix is None:
             if max_supercell is None:
-                raise ValueError("If scale_matrix is not provided, please provide max_supercell!")
-
+                raise ValueError("If scale_matrix is not provided in input, please provide max_supercell!")
             scale_matrix = self.get_smallest_supercell(qpoint, max_supercell=max_supercell)
 
         scale_matrix = np.array(scale_matrix, np.int16)
         if scale_matrix.shape != (3, 3):
             scale_matrix = np.array(scale_matrix * np.eye(3), np.int16)
-        print("scale_matrix:", scale_matrix)
+        #print("scale_matrix:\n", scale_matrix)
 
         old_lattice = self._lattice
         new_lattice = Lattice(np.dot(scale_matrix, old_lattice.matrix))
 
         tvects = self.get_trans_vect(scale_matrix)
-        print("tvects", tvects)
+        #print("tvects\n", tvects)
 
         if frac_coords:
-            displ = np.array((old_lattice.get_cartesian_coords(d) for d in displ))
+            displ = np.array([old_lattice.get_cartesian_coords(d) for d in displ])
         else:
             displ = np.array(displ)
-        # from here displ are in cartesian coordinates
 
+        # from here on, displ is in cartesian coordinates.
         displ = eta * displ / np.linalg.norm(displ, axis=1).max()
 
         new_displ = np.zeros(3, dtype=float)
         new_sites = []
         displ_list = []
-        for at, site in enumerate(self):
+        for iat, site in enumerate(self):
             for t in tvects:
-                new_displ[:] = np.real(np.exp(2*1j*np.pi*(np.dot(qpoint,t)))*displ[at,:])
+                new_displ[:] = np.real(np.exp(2*1j*np.pi * (np.dot(qpoint, t))) * displ[iat,:])
                 displ_list.append(list(new_displ))
 
                 coords = site.coords + old_lattice.get_cartesian_coords(t) + new_displ
@@ -2699,7 +2738,7 @@ def structure2siesta(structure: Structure, verbose=0) -> str:
     lines = []
     app = lines.append
     app("NumberOfAtoms %d" % len(structure))
-    app("NumberOfSpecies %d" % structure.ntypesp)
+    app("NumberOfSpecies %d" % structure.n_elems)
 
     if verbose:
         app("# The species number followed by the atomic number, and then by the desired label")

abipy/core/testing.py

@@ -39,7 +39,7 @@ def cmp_version(this: str, other: str, op: str = ">=") -> bool:
     Compare two version strings with the given operator ``op``
     >>> assert cmp_version("1.1.1", "1.1.0") and not cmp_version("1.1.1", "1.1.0", op="==")
     """
-    from pkg_resources import parse_version
+    from packaging.version import parse as parse_version
     from monty.operator import operator_from_str
     op = operator_from_str(op)
     return op(parse_version(this), parse_version(other))

abipy/core/tests/test_structure.py

@@ -397,7 +397,7 @@ def test_frozen_phonon_methods(self):
         f2p_data = structure.frozen_2phonon(qpoint, 0.05 * displ, 0.02*displ2, eta=0.5, frac_coords=False,
                                             max_supercell=mx_sc, scale_matrix=scale_matrix)
 
-        d_tot = 0.05*displ+0.02*displ2
+        d_tot = 0.05 *displ + 0.02 * displ2
         max_displ = np.linalg.norm(d_tot, axis=1).max()
         self.assert_almost_equal(f2p_data.structure[0].coords,
                                     structure[0].coords + 0.5*d_tot[0]/max_displ)