From 187d8752eebcc4158090ef7b3058b8d0f37e7466 Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Sun, 21 Jan 2024 09:23:48 +0100 Subject: [PATCH 1/5] first try at implementing bond-order calculations This fixes #507 by adding the Wiberg and Mayer implementations. I am not fully sure the Mayer is correct, since the expressions are PS_abPS_ba, but due to symmetry, I would expect this to be PS_ab ** 2. fixed bond+antibond method, factor 0.5 Signed-off-by: Nick Papior --- CHANGELOG.md | 1 + src/sisl/_core/sparse.py | 42 ++++-- src/sisl/io/siesta/_help.py | 1 + src/sisl/physics/densitymatrix.py | 125 +++++++++++++++++- src/sisl/physics/tests/test_density_matrix.py | 9 ++ 5 files changed, 166 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index f4597b0970..d52d47ab9a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ we hit release version 1.0.0. ### Added - added an efficient neighbor finder, #393 - enabled reading DFTB+ output Hamiltonian and overlap matrices, #579 +- `bond_order` for `DensityMatrix` objects, #507 - better error messages when users request quantities not calculated by Siesta/TBtrans - functional programming of the basic sisl classes Now many of the `Geometry|Lattice|Grid.* manipulation routines which diff --git a/src/sisl/_core/sparse.py b/src/sisl/_core/sparse.py index cc00049b53..fffd53a50f 100644 --- a/src/sisl/_core/sparse.py +++ b/src/sisl/_core/sparse.py @@ -57,24 +57,44 @@ __all__ = ["SparseCSR", "ispmatrix", "ispmatrixd"] -def _rows_and_cols(csr, rows=None): - """Retrieve the sparse patterns rows and columns from the sparse matrix +def _to_coo(csr, data: bool = True, rows=None): + """Retrieve the CSR information in a sanitized manner - Possibly only retrieving it for a subset of rows. + I.e. as COO format with data stripped, or only rows and columns + + Parameters + ---------- + csr: SparseCSR + matrix to sanitize + data: bool + whether the data should also be returned sanitized + rows: + only return for a subset of rowsj + + Returns + ------- + rows, cols : always + matrix_data : when `data` is True """ ptr = csr.ptr ncol = csr.ncol col = csr.col + D = csr._D + if rows is None: - cols = col[array_arange(ptr[:-1], n=ncol, dtype=int32)] - idx = (ncol > 0).nonzero()[0] - rows = repeat(idx.astype(int32, copy=False), ncol[idx]) + idx = array_arange(ptr[:-1], n=ncol, dtype=int32) else: rows = csr._sanitize(rows).ravel() ncol = ncol[rows] - cols = col[array_arange(ptr[rows], n=ncol, dtype=int32)] - idx = (ncol > 0).nonzero()[0] - rows = repeat(rows[idx].astype(int32, copy=False), ncol[idx]) + idx = array_arange(ptr[rows], n=ncol, dtype=int32) + if data: + D = D[idx] + cols = col[idx] + idx = (ncol > 0).nonzero()[0] + rows = repeat(idx.astype(int32, copy=False), ncol[idx]) + + if data: + return rows, cols, D return rows, cols @@ -342,7 +362,7 @@ def diagonal(self): # get the diagonal components diag = np.zeros([self.shape[0], self.shape[2]], dtype=self.dtype) - rows, cols = _rows_and_cols(self) + rows, cols = _to_coo(self, data=False) # Now retrieve rows and cols idx = array_arange(self.ptr[:-1], n=self.ncol, dtype=int32) @@ -2011,7 +2031,7 @@ def rowslice(r): accessors = mat.dim, mat.col, mat._D, rowslice # check whether they are actually sorted if mat.finalized: - rows, cols = _rows_and_cols(mat) + rows, cols = _to_coo(mat, data=False) rows, cols = np.diff(rows), np.diff(cols) issorted = np.all(cols[rows == 0] > 0) else: diff --git a/src/sisl/io/siesta/_help.py b/src/sisl/io/siesta/_help.py index 3884831a94..64aa5a89e7 100644 --- a/src/sisl/io/siesta/_help.py +++ b/src/sisl/io/siesta/_help.py @@ -7,6 +7,7 @@ from sisl import SislError from sisl._core.sparse import _rows_and_cols from sisl.messages import warn +from sisl.sparse import _to_coo try: from . import _siesta diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py index d56f71474f..29cd90bcf1 100644 --- a/src/sisl/physics/densitymatrix.py +++ b/src/sisl/physics/densitymatrix.py @@ -6,7 +6,7 @@ import numpy as np from numpy import add, dot, logical_and, repeat, subtract, unique -from scipy.sparse import csr_matrix +from scipy.sparse import coo_matrix, csr_matrix from scipy.sparse import hstack as ss_hstack from scipy.sparse import tril, triu @@ -19,6 +19,8 @@ from sisl._internal import set_module from sisl._math_small import xyz_to_spherical_cos_phi from sisl.messages import progressbar, warn +from sisl.sparse import SparseCSR, _ncol_to_indptr, _to_coo +from sisl.sparse_geometry import SparseAtom, SparseOrbital from .sparse import SparseOrbitalBZSpin from .spin import Spin @@ -26,6 +28,30 @@ __all__ = ["DensityMatrix"] +def _get_density(DM, what="sum"): + DM = DM.T + if what == "sum": + if DM.shape[0] in (2, 4, 8): + return DM[0] + DM[1] + return DM[0] + if what == "spin": + m = np.empty([3, DM.shape[1]], dtype=DM.dtype) + if DM.shape[0] == 8: + m[0] = DM[2] + DM[6] + m[1] = -DM[3] + DM[7] + m[2] = DM[0] - DM[1] + elif DM.shape[0] == 4: + m[0] = 2 * DM[2] + m[1] = -2 * DM[3] + m[2] = DM[0] - DM[1] + elif DM.shape[0] == 2: + m[:2, :] = 0.0 + m[2] = DM[0] - DM[1] + elif DM.shape[0] == 1: + m[...] = 0.0 + return m + + class _densitymatrix(SparseOrbitalBZSpin): def spin_rotate(self, angles, rad=False): r"""Rotates spin-boxes by fixed angles around the :math:`x`, :math:`y` and :math:`z` axis, respectively. @@ -397,6 +423,103 @@ def _convert(M): f"{self.__class__.__name__}.mulliken only allows projection [orbital, atom]" ) + def bond_order(self, method: str = "mayer"): + r"""Bond-order calculation using Mayer, Wiberg or other methods + + For ``method='wiberg'``, the bond-order is calculated as: + + .. math:: + B_{\alpha\beta}^{\mathrm{Wiberg}} = \sum_{\nu\in \alpha}\sum_{\mu\in \beta} D_{\nu\mu}^2 + + For ``method='mayer'``, the bond-order is calculated as: + + .. math:: + B_{\alpha\beta}^{\mathrm{Mayer}} = \sum_{\nu\in \alpha}\sum_{\mu\in \beta} D_{\nu\mu}S_{\nu\mu}D_{\mu\nu}S_{\mu\nu} + + + For ``method='bond+anti'``, the bond-order is calculated as: + + .. math:: + B_{\alpha\beta}^{\mathrm{b}+\mathrm{ab}} = \sum_{\nu\in \alpha}\sum_{\mu\in \beta} D_{\nu\mu}S_{\nu\mu} / 2 + + For all options one can do the bond-order calculation for the + spin components. Albeit, their meaning may be more doubtful. + Simply add ``':spin'`` to the `method` argument, and the returned + quantity will be spin-resolved. + + + Note + ---- + Wiberg and Mayer will be equivalent for orthogonal basis sets. + + It is unclear what the spin-density bond-order really means. + + Parameters + ---------- + method : {mayer, wiberg, bond+anti}[:spin] + which method to calculate the bond-order with + + Returns + ------- + SparseAtom : with the bond-order between any two atoms, in a supercell matrix. + """ + method = method.lower() + + # split method to retrieve options + m, *opts = method.split(":") + + # only extract the summed density + what = "sum" + if "spin" in opts: + # do this for each spin x, y, z + what = "spin" + del opts[opts.index("spin")] + + # Check that there are no un-used options + if opts: + raise ValueError( + f"{self.__class__.__name__}.bond_order got non-valid options {opts}" + ) + + # get all rows and columns + rows, cols, DM = _to_coo(self._csr) + + # Add factor S, if needed + if not self.orthogonal: + if method != "wiberg": + DM = DM[:, :-1] * DM[:, -1].reshape(-1, 1) + + # Convert to requested matrix form + DM = _get_density(DM, what) + if m in ("wiberg", "mayer"): + # calculate the bond-order matrix + DM **= 2 + elif m == "bond+anti": + # sum of COOP for each atom + DM /= 2 + else: + raise ValueError( + f"{self.__class__.__name__}.bond_order got non-valid method {method}" + ) + + # create the sparseatom object + geom = self.geometry + rows = geom.o2a(rows) + cols = geom.o2a(cols) + + shape = (geom.na, geom.na_s) + if DM.ndim == 2: + data = [] + for d in DM: + d = coo_matrix((d, (rows, cols)), shape=shape).tocsr() + d.sum_duplicates() + data.append(d) + else: + data = coo_matrix((DM, (rows, cols)), shape=shape).tocsr() + data.sum_duplicates() + + return SparseAtom.fromsp(geom, data) + def density(self, grid, spinor=None, tol=1e-7, eta=None): r"""Expand the density matrix to the charge density on a grid diff --git a/src/sisl/physics/tests/test_density_matrix.py b/src/sisl/physics/tests/test_density_matrix.py index 2bbd0bd67f..1850f83641 100644 --- a/src/sisl/physics/tests/test_density_matrix.py +++ b/src/sisl/physics/tests/test_density_matrix.py @@ -13,6 +13,7 @@ Geometry, Grid, Lattice, + SparseAtom, SphericalOrbital, Spin, ) @@ -165,6 +166,14 @@ def test_mulliken_polarized(self): assert m[0].sum() == pytest.approx(3) assert m[1].sum() == pytest.approx(1) + @pytest.mark.parametrize("method", ["wiberg", "mayer", "bond+anti"]) + @pytest.mark.parametrize("option", ["", ":spin"]) + def test_bond_order(self, setup, method, option): + D = setup.D.copy() + D.construct(setup.func) + BO = D.bond_order(method + option) + assert isinstance(BO, SparseAtom) + def test_rho1(self, setup): D = setup.D.copy() D.construct(setup.func) From 54e2cf69f2ed1eb2c4b86813d423b6c3df4603ad Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Tue, 27 Feb 2024 13:44:08 +0100 Subject: [PATCH 2/5] fixed bond-order calculations for wiberg and mayer Still needs a bit refinement, but it should be there now. Signed-off-by: Nick Papior --- src/sisl/io/siesta/_help.py | 2 - src/sisl/physics/densitymatrix.py | 154 +++++++++++++----- src/sisl/physics/tests/test_density_matrix.py | 3 +- 3 files changed, 118 insertions(+), 41 deletions(-) diff --git a/src/sisl/io/siesta/_help.py b/src/sisl/io/siesta/_help.py index 64aa5a89e7..4e340cba12 100644 --- a/src/sisl/io/siesta/_help.py +++ b/src/sisl/io/siesta/_help.py @@ -5,9 +5,7 @@ import sisl._array as _a from sisl import SislError -from sisl._core.sparse import _rows_and_cols from sisl.messages import warn -from sisl.sparse import _to_coo try: from . import _siesta diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py index 29cd90bcf1..40f8f49301 100644 --- a/src/sisl/physics/densitymatrix.py +++ b/src/sisl/physics/densitymatrix.py @@ -13,14 +13,12 @@ import sisl._array as _a from sisl import BoundaryCondition as BC from sisl import Geometry, Lattice -from sisl._core.sparse import SparseCSR, _ncol_to_indptr -from sisl._core.sparse_geometry import SparseOrbital +from sisl._core.sparse import SparseCSR, _ncol_to_indptr, _to_coo +from sisl._core.sparse_geometry import SparseAtom, SparseOrbital from sisl._indices import indices_fabs_le, indices_le from sisl._internal import set_module from sisl._math_small import xyz_to_spherical_cos_phi from sisl.messages import progressbar, warn -from sisl.sparse import SparseCSR, _ncol_to_indptr, _to_coo -from sisl.sparse_geometry import SparseAtom, SparseOrbital from .sparse import SparseOrbitalBZSpin from .spin import Spin @@ -28,26 +26,30 @@ __all__ = ["DensityMatrix"] -def _get_density(DM, what="sum"): +def _get_density(DM, orthogonal, what="sum"): DM = DM.T + if orthogonal: + off = 0 + else: + off = 1 if what == "sum": - if DM.shape[0] in (2, 4, 8): + if DM.shape[0] in (2 + off, 4 + off, 8 + off): return DM[0] + DM[1] return DM[0] if what == "spin": m = np.empty([3, DM.shape[1]], dtype=DM.dtype) - if DM.shape[0] == 8: + if DM.shape[0] == 8 + off: m[0] = DM[2] + DM[6] m[1] = -DM[3] + DM[7] m[2] = DM[0] - DM[1] - elif DM.shape[0] == 4: + elif DM.shape[0] == 4 + off: m[0] = 2 * DM[2] m[1] = -2 * DM[3] m[2] = DM[0] - DM[1] - elif DM.shape[0] == 2: + elif DM.shape[0] == 2 + off: m[:2, :] = 0.0 m[2] = DM[0] - DM[1] - elif DM.shape[0] == 1: + elif DM.shape[0] == 1 + off: m[...] = 0.0 return m @@ -482,44 +484,120 @@ def bond_order(self, method: str = "mayer"): ) # get all rows and columns + geom = self.geometry rows, cols, DM = _to_coo(self._csr) - # Add factor S, if needed - if not self.orthogonal: - if method != "wiberg": - DM = DM[:, :-1] * DM[:, -1].reshape(-1, 1) + # Define a matrix-matrix multiplication + def mm(A, B): + n = A.shape[0] + latt = self.geometry.lattice + sc_off = latt.sc_off + + # A matrix product in a supercell is a bit tricky + # since the off-diagonal elements are formed with + # respect to the supercell offsets from the diagonal + # compoent + + # A = [[ sc1-sc1, sc2-sc1, sc3-sc1, ... + # sc1-sc2, sc2-sc2, sc3-sc2, ... + # sc1-sc3, sc2-sc3, sc3-sc3, ... + + # so each column has a *principal* supercell + # which is used to calculate the offset of each + # other component. + # Now for the LHS in a MM, we have A[0, :] + # which is only wrt. the 0,0 component. + # In sisl this is forced to be the supercell 0,0. + # Hence everything in that row requires no special + # handling. Yet all others do. + + res = [] + for i_s in range(latt.n_s): + + # Calculate the 0,i_s column of the MM + # This is equal to: + # A[0, :] @ B[:, i_s] + # Calculate the offset for the B column + sc_offj = sc_off[i_s] - sc_off + + # initialize the result array + # Not strictly needed, but enforces that the + # data always contains a csr_matrix + r = csr_matrix((n, n), dtype=A.dtype) + + # get current supercell information + for i in range(latt.n_s): + # i == LHS matrix + # j == RHS matrix + try: + sc = sc_offj[i] + # if the supercell index does not exist, it means + # the matrix is 0. Hence we just neglect that contribution. + j = latt.sc_index(sc) + r = r + A[:, i * n : (i + 1) * n] @ B[:, j * n : (j + 1) * n] + except: + continue + + res.append(r) + + # Re-create a matrix where each block is joined into a + # big matrix, hstack == columnwise stacking. + return ss_hstack(res) - # Convert to requested matrix form - DM = _get_density(DM, what) if m in ("wiberg", "mayer"): - # calculate the bond-order matrix - DM **= 2 + if self.orthogonal: + S = np.zeros(rows.size, dtype=DM.dtype) + S[rows == cols] = 1.0 + else: + S = DM[:, -1] + + # Convert to requested matrix form + D = _get_density(DM, self.orthogonal, what) + + def get_BO(geom, D, S, rows, cols): + # TODO, for each spin-component, do this: + DM = self.fromsp( + geom, + coo_matrix((D, (rows, cols)), shape=self.shape[:2]), + S=coo_matrix((S, (rows, cols)), shape=self.shape[:2]), + ) + + D = DM.Dk(format="sc:csr") + if m == "mayer": + S = DM.Sk(format="sc:csr") + BO = mm(D, S).multiply(mm(S, D)) + else: + # we should probably do something else + # the Wiberg is sum_A sum_B D_ab + # hence a double summation would be needed. + # Currently we just multiply by 2 + # Probably wrong, since we need the off-diagonal + BO = mm(D, D) * 2 + + BO = BO.tocoo() + + # Now re-create the sparse-atom component. + rows = geom.o2a(BO.row) + cols = geom.o2a(BO.col) + shape = (geom.na, geom.na_s) + BO = coo_matrix((BO.data, (rows, cols)), shape=shape).tocsr() + BO.sum_duplicates() + return BO + + if D.ndim == 2: + BO = [get_BO(geom, d, S, rows, cols) for d in D] + else: + BO = get_BO(geom, D, S, rows, cols) + + return SparseAtom.fromsp(geom, BO) + elif m == "bond+anti": - # sum of COOP for each atom - DM /= 2 + raise NotImplementedError else: raise ValueError( f"{self.__class__.__name__}.bond_order got non-valid method {method}" ) - # create the sparseatom object - geom = self.geometry - rows = geom.o2a(rows) - cols = geom.o2a(cols) - - shape = (geom.na, geom.na_s) - if DM.ndim == 2: - data = [] - for d in DM: - d = coo_matrix((d, (rows, cols)), shape=shape).tocsr() - d.sum_duplicates() - data.append(d) - else: - data = coo_matrix((DM, (rows, cols)), shape=shape).tocsr() - data.sum_duplicates() - - return SparseAtom.fromsp(geom, data) - def density(self, grid, spinor=None, tol=1e-7, eta=None): r"""Expand the density matrix to the charge density on a grid diff --git a/src/sisl/physics/tests/test_density_matrix.py b/src/sisl/physics/tests/test_density_matrix.py index 1850f83641..06ef8bfa79 100644 --- a/src/sisl/physics/tests/test_density_matrix.py +++ b/src/sisl/physics/tests/test_density_matrix.py @@ -87,6 +87,7 @@ def func(D, ia, atoms, atoms_xyz): @pytest.mark.physics @pytest.mark.density_matrix +@pytest.mark.densitymatrix class TestDensityMatrix: def test_objects(self, setup): assert len(setup.D.xyz) == 2 @@ -166,7 +167,7 @@ def test_mulliken_polarized(self): assert m[0].sum() == pytest.approx(3) assert m[1].sum() == pytest.approx(1) - @pytest.mark.parametrize("method", ["wiberg", "mayer", "bond+anti"]) + @pytest.mark.parametrize("method", ["wiberg", "mayer"]) @pytest.mark.parametrize("option", ["", ":spin"]) def test_bond_order(self, setup, method, option): D = setup.D.copy() From ee5fd312bc8c5aa32c86d3a84c37a84ebcf12eb1 Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Wed, 28 Feb 2024 12:14:51 +0100 Subject: [PATCH 3/5] speeded up large system calculations Also corrected Wiberg (which gives bad results) Signed-off-by: Nick Papior --- src/sisl/physics/densitymatrix.py | 82 +++++++++++++++++-------------- 1 file changed, 45 insertions(+), 37 deletions(-) diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py index 40f8f49301..250a8230a8 100644 --- a/src/sisl/physics/densitymatrix.py +++ b/src/sisl/physics/densitymatrix.py @@ -431,12 +431,12 @@ def bond_order(self, method: str = "mayer"): For ``method='wiberg'``, the bond-order is calculated as: .. math:: - B_{\alpha\beta}^{\mathrm{Wiberg}} = \sum_{\nu\in \alpha}\sum_{\mu\in \beta} D_{\nu\mu}^2 + B_{\alpha\beta}^{\mathrm{Wiberg}} = \sum_{\nu\in\alpha}\sum_{\mu\in\beta} D_{\nu\mu}^2 For ``method='mayer'``, the bond-order is calculated as: .. math:: - B_{\alpha\beta}^{\mathrm{Mayer}} = \sum_{\nu\in \alpha}\sum_{\mu\in \beta} D_{\nu\mu}S_{\nu\mu}D_{\mu\nu}S_{\mu\nu} + B_{\alpha\beta}^{\mathrm{Mayer}} = \sum_{\nu\in\alpha}\sum_{\mu\in\beta} (DS)_{\nu\mu}(DS)_{\mu\nu} For ``method='bond+anti'``, the bond-order is calculated as: @@ -452,8 +452,6 @@ def bond_order(self, method: str = "mayer"): Note ---- - Wiberg and Mayer will be equivalent for orthogonal basis sets. - It is unclear what the spin-density bond-order really means. Parameters @@ -487,12 +485,20 @@ def bond_order(self, method: str = "mayer"): geom = self.geometry rows, cols, DM = _to_coo(self._csr) + # Convert to requested matrix form + D = _get_density(DM, self.orthogonal, what) + # Define a matrix-matrix multiplication def mm(A, B): n = A.shape[0] latt = self.geometry.lattice sc_off = latt.sc_off + # we will extract sub-matrices n_s ** 2 times. + # Extracting once should be fine (and ok) + Al = [A[:, i * n : (i + 1) * n] for i in range(latt.n_s)] + Bl = [B[:, i * n : (i + 1) * n] for i in range(latt.n_s)] + # A matrix product in a supercell is a bit tricky # since the off-diagonal elements are formed with # respect to the supercell offsets from the diagonal @@ -534,55 +540,57 @@ def mm(A, B): # if the supercell index does not exist, it means # the matrix is 0. Hence we just neglect that contribution. j = latt.sc_index(sc) - r = r + A[:, i * n : (i + 1) * n] @ B[:, j * n : (j + 1) * n] + r = r + Al[i] @ Bl[j] except: continue res.append(r) + # Clean-up... + del Al, Bl + # Re-create a matrix where each block is joined into a # big matrix, hstack == columnwise stacking. return ss_hstack(res) - if m in ("wiberg", "mayer"): + def sparseorb2sparseatom(geom, M): + # Ensure we have in COO-rdinate format + M = M.tocoo() + + # Now re-create the sparse-atom component. + rows = geom.o2a(M.row) + cols = geom.o2a(M.col) + shape = (geom.na, geom.na_s) + M = coo_matrix((M.data, (rows, cols)), shape=shape).tocsr() + M.sum_duplicates() + return M + + if m == "wiberg": + + def get_BO(geom, D, rows, cols): + # square of each element + BO = coo_matrix((D * D, (rows, cols)), shape=self.shape[:2]) + return sparseorb2sparseatom(geom, BO) + + if D.ndim == 2: + BO = [get_BO(geom, d, rows, cols) for d in D] + else: + BO = get_BO(geom, D, rows, cols) + + return SparseAtom.fromsp(geom, BO) + + elif m == "mayer": if self.orthogonal: S = np.zeros(rows.size, dtype=DM.dtype) S[rows == cols] = 1.0 else: S = DM[:, -1] - # Convert to requested matrix form - D = _get_density(DM, self.orthogonal, what) - def get_BO(geom, D, S, rows, cols): - # TODO, for each spin-component, do this: - DM = self.fromsp( - geom, - coo_matrix((D, (rows, cols)), shape=self.shape[:2]), - S=coo_matrix((S, (rows, cols)), shape=self.shape[:2]), - ) - - D = DM.Dk(format="sc:csr") - if m == "mayer": - S = DM.Sk(format="sc:csr") - BO = mm(D, S).multiply(mm(S, D)) - else: - # we should probably do something else - # the Wiberg is sum_A sum_B D_ab - # hence a double summation would be needed. - # Currently we just multiply by 2 - # Probably wrong, since we need the off-diagonal - BO = mm(D, D) * 2 - - BO = BO.tocoo() - - # Now re-create the sparse-atom component. - rows = geom.o2a(BO.row) - cols = geom.o2a(BO.col) - shape = (geom.na, geom.na_s) - BO = coo_matrix((BO.data, (rows, cols)), shape=shape).tocsr() - BO.sum_duplicates() - return BO + D = coo_matrix((D, (rows, cols)), shape=self.shape[:2]).tocsr() + S = coo_matrix((S, (rows, cols)), shape=self.shape[:2]).tocsr() + BO = mm(D, S).multiply(mm(S, D)) + return sparseorb2sparseatom(geom, BO) if D.ndim == 2: BO = [get_BO(geom, d, S, rows, cols) for d in D] From 40907daaa0f5142e555fd90687afb9ba213c675e Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Wed, 28 Feb 2024 21:30:05 +0100 Subject: [PATCH 4/5] added Mulliken bond-order While the mulliken bond-order can easily be calculated from the mulliken output, this will make it simpler. Signed-off-by: Nick Papior --- src/sisl/physics/densitymatrix.py | 62 ++++++++++++++++--------------- 1 file changed, 33 insertions(+), 29 deletions(-) diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py index 250a8230a8..5640dab39c 100644 --- a/src/sisl/physics/densitymatrix.py +++ b/src/sisl/physics/densitymatrix.py @@ -426,7 +426,7 @@ def _convert(M): ) def bond_order(self, method: str = "mayer"): - r"""Bond-order calculation using Mayer, Wiberg or other methods + r"""Bond-order calculation using various methods For ``method='wiberg'``, the bond-order is calculated as: @@ -438,17 +438,16 @@ def bond_order(self, method: str = "mayer"): .. math:: B_{\alpha\beta}^{\mathrm{Mayer}} = \sum_{\nu\in\alpha}\sum_{\mu\in\beta} (DS)_{\nu\mu}(DS)_{\mu\nu} - - For ``method='bond+anti'``, the bond-order is calculated as: + For ``method='mulliken'``, the bond-order is calculated as: .. math:: - B_{\alpha\beta}^{\mathrm{b}+\mathrm{ab}} = \sum_{\nu\in \alpha}\sum_{\mu\in \beta} D_{\nu\mu}S_{\nu\mu} / 2 + B_{\alpha\beta}^{\mathrm{Mulliken}} = 2\sum_{\nu\in\alpha}\sum_{\mu\in\beta} D_{\nu\mu}S_{\nu\mu} For all options one can do the bond-order calculation for the spin components. Albeit, their meaning may be more doubtful. Simply add ``':spin'`` to the `method` argument, and the returned - quantity will be spin-resolved. - + quantity will be spin-resolved with :math:`x`, :math:`y` and :math:`z` + components. Note ---- @@ -456,7 +455,7 @@ def bond_order(self, method: str = "mayer"): Parameters ---------- - method : {mayer, wiberg, bond+anti}[:spin] + method : {mayer, wiberg, mulliken}[:spin] which method to calculate the bond-order with Returns @@ -532,11 +531,10 @@ def mm(A, B): r = csr_matrix((n, n), dtype=A.dtype) # get current supercell information - for i in range(latt.n_s): + for i, sc in enumerate(sc_offj): # i == LHS matrix # j == RHS matrix try: - sc = sc_offj[i] # if the supercell index does not exist, it means # the matrix is 0. Hence we just neglect that contribution. j = latt.sc_index(sc) @@ -565,26 +563,18 @@ def sparseorb2sparseatom(geom, M): M.sum_duplicates() return M + S = False + if m == "wiberg": - def get_BO(geom, D, rows, cols): + def get_BO(geom, D, S, rows, cols): # square of each element BO = coo_matrix((D * D, (rows, cols)), shape=self.shape[:2]) return sparseorb2sparseatom(geom, BO) - if D.ndim == 2: - BO = [get_BO(geom, d, rows, cols) for d in D] - else: - BO = get_BO(geom, D, rows, cols) - - return SparseAtom.fromsp(geom, BO) - elif m == "mayer": - if self.orthogonal: - S = np.zeros(rows.size, dtype=DM.dtype) - S[rows == cols] = 1.0 - else: - S = DM[:, -1] + + S = True def get_BO(geom, D, S, rows, cols): D = coo_matrix((D, (rows, cols)), shape=self.shape[:2]).tocsr() @@ -592,20 +582,34 @@ def get_BO(geom, D, S, rows, cols): BO = mm(D, S).multiply(mm(S, D)) return sparseorb2sparseatom(geom, BO) - if D.ndim == 2: - BO = [get_BO(geom, d, S, rows, cols) for d in D] - else: - BO = get_BO(geom, D, S, rows, cols) + elif m == "mulliken": + + S = True - return SparseAtom.fromsp(geom, BO) + def get_BO(geom, D, S, rows, cols): + # Got the factor 2 from Multiwfn + BO = coo_matrix((D * S * 2, (rows, cols)), shape=self.shape[:2]).tocsr() + return sparseorb2sparseatom(geom, BO) - elif m == "bond+anti": - raise NotImplementedError else: raise ValueError( f"{self.__class__.__name__}.bond_order got non-valid method {method}" ) + if S: + if self.orthogonal: + S = np.zeros(rows.size, dtype=DM.dtype) + S[rows == cols] = 1.0 + else: + S = DM[:, -1] + + if D.ndim == 2: + BO = [get_BO(geom, d, S, rows, cols) for d in D] + else: + BO = get_BO(geom, D, S, rows, cols) + + return SparseAtom.fromsp(geom, BO) + def density(self, grid, spinor=None, tol=1e-7, eta=None): r"""Expand the density matrix to the charge density on a grid From b65bbe851c468ad3e51ce65aaff32822c8a7cfd2 Mon Sep 17 00:00:00 2001 From: Nick Papior Date: Thu, 29 Feb 2024 09:28:31 +0100 Subject: [PATCH 5/5] added orbital projection to bond-order While it is a bit weird to have bond_order return orbital-order, I think it is not necessary to create a new method to retrieve the bond-order. Better to keep it simple. Signed-off-by: Nick Papior --- src/sisl/physics/densitymatrix.py | 76 ++++++++++++++----- src/sisl/physics/tests/test_density_matrix.py | 13 +++- 2 files changed, 65 insertions(+), 24 deletions(-) diff --git a/src/sisl/physics/densitymatrix.py b/src/sisl/physics/densitymatrix.py index 5640dab39c..0f190c5456 100644 --- a/src/sisl/physics/densitymatrix.py +++ b/src/sisl/physics/densitymatrix.py @@ -425,23 +425,28 @@ def _convert(M): f"{self.__class__.__name__}.mulliken only allows projection [orbital, atom]" ) - def bond_order(self, method: str = "mayer"): + def bond_order(self, method: str = "mayer", projection: str = "atom"): r"""Bond-order calculation using various methods For ``method='wiberg'``, the bond-order is calculated as: .. math:: - B_{\alpha\beta}^{\mathrm{Wiberg}} = \sum_{\nu\in\alpha}\sum_{\mu\in\beta} D_{\nu\mu}^2 + B_{\nu\mu}^{\mathrm{Wiberg}} &= D_{\nu\mu}^2 + \\ + B_{\alpha\beta}^{\mathrm{Wiberg}} &= \sum_{\nu\in\alpha}\sum_{\mu\in\beta} B_{\nu\mu} For ``method='mayer'``, the bond-order is calculated as: .. math:: - B_{\alpha\beta}^{\mathrm{Mayer}} = \sum_{\nu\in\alpha}\sum_{\mu\in\beta} (DS)_{\nu\mu}(DS)_{\mu\nu} - - For ``method='mulliken'``, the bond-order is calculated as: - + B_{\nu\mu}^{\mathrm{Mayer}} &= (DS)_{\nu\mu}(DS)_{\mu\nu} + \\ + B_{\alpha\beta}^{\mathrm{Mayer}} &= \sum_{\nu\in\alpha}\sum_{\mu\in\beta} B_{\nu\mu} .. math:: - B_{\alpha\beta}^{\mathrm{Mulliken}} = 2\sum_{\nu\in\alpha}\sum_{\mu\in\beta} D_{\nu\mu}S_{\nu\mu} + B_{\nu\mu}^{\mathrm{Mulliken}} &= 2D_{\nu\mu}S{\nu\mu} + \\ + B_{\alpha\beta}^{\mathrm{Mulliken}} &= \sum_{\nu\in\alpha}\sum_{\mu\in\beta} B__{\nu\mu} + + The Mulliken bond-order is closely related to the COOP interpretation. For all options one can do the bond-order calculation for the spin components. Albeit, their meaning may be more doubtful. @@ -458,9 +463,17 @@ def bond_order(self, method: str = "mayer"): method : {mayer, wiberg, mulliken}[:spin] which method to calculate the bond-order with + projection : {atom, orbital} + whether the returned matrix is in orbital form, or in atom form. + If orbital is used, then the above formulas will be changed + + Returns ------- SparseAtom : with the bond-order between any two atoms, in a supercell matrix. + Returned only if projection is atom. + SparseOrbital : with the bond-order between any two orbitals, in a supercell matrix. + Returned only if projection is orbital. """ method = method.lower() @@ -551,17 +564,38 @@ def mm(A, B): # big matrix, hstack == columnwise stacking. return ss_hstack(res) - def sparseorb2sparseatom(geom, M): - # Ensure we have in COO-rdinate format - M = M.tocoo() + projection = projection.lower() + + if projection.startswith("atom"): # allows atoms + + out_cls = SparseAtom + + def sparse2sparse(geom, M): + + # Ensure we have in COO-rdinate format + M = M.tocoo() - # Now re-create the sparse-atom component. - rows = geom.o2a(M.row) - cols = geom.o2a(M.col) - shape = (geom.na, geom.na_s) - M = coo_matrix((M.data, (rows, cols)), shape=shape).tocsr() - M.sum_duplicates() - return M + # Now re-create the sparse-atom component. + rows = geom.o2a(M.row) + cols = geom.o2a(M.col) + shape = (geom.na, geom.na_s) + M = coo_matrix((M.data, (rows, cols)), shape=shape).tocsr() + M.sum_duplicates() + return M + + elif projection.startswith("orbital"): # allows orbitals + + out_cls = SparseOrbital + + def sparse2sparse(geom, M): + M = M.tocsr() + M.sum_duplicates() + return M + + else: + raise ValueError( + f"{self.__class__.__name__}.bond_order got unexpected keyword projection" + ) S = False @@ -570,7 +604,7 @@ def sparseorb2sparseatom(geom, M): def get_BO(geom, D, S, rows, cols): # square of each element BO = coo_matrix((D * D, (rows, cols)), shape=self.shape[:2]) - return sparseorb2sparseatom(geom, BO) + return sparse2sparse(geom, BO) elif m == "mayer": @@ -580,7 +614,7 @@ def get_BO(geom, D, S, rows, cols): D = coo_matrix((D, (rows, cols)), shape=self.shape[:2]).tocsr() S = coo_matrix((S, (rows, cols)), shape=self.shape[:2]).tocsr() BO = mm(D, S).multiply(mm(S, D)) - return sparseorb2sparseatom(geom, BO) + return sparse2sparse(geom, BO) elif m == "mulliken": @@ -589,7 +623,7 @@ def get_BO(geom, D, S, rows, cols): def get_BO(geom, D, S, rows, cols): # Got the factor 2 from Multiwfn BO = coo_matrix((D * S * 2, (rows, cols)), shape=self.shape[:2]).tocsr() - return sparseorb2sparseatom(geom, BO) + return sparse2sparse(geom, BO) else: raise ValueError( @@ -608,7 +642,7 @@ def get_BO(geom, D, S, rows, cols): else: BO = get_BO(geom, D, S, rows, cols) - return SparseAtom.fromsp(geom, BO) + return out_cls.fromsp(geom, BO) def density(self, grid, spinor=None, tol=1e-7, eta=None): r"""Expand the density matrix to the charge density on a grid diff --git a/src/sisl/physics/tests/test_density_matrix.py b/src/sisl/physics/tests/test_density_matrix.py index 06ef8bfa79..a5a4a19b0c 100644 --- a/src/sisl/physics/tests/test_density_matrix.py +++ b/src/sisl/physics/tests/test_density_matrix.py @@ -14,6 +14,7 @@ Grid, Lattice, SparseAtom, + SparseOrbital, SphericalOrbital, Spin, ) @@ -169,11 +170,17 @@ def test_mulliken_polarized(self): @pytest.mark.parametrize("method", ["wiberg", "mayer"]) @pytest.mark.parametrize("option", ["", ":spin"]) - def test_bond_order(self, setup, method, option): + @pytest.mark.parametrize("projection", ["atom", "orbitals"]) + def test_bond_order(self, setup, method, option, projection): D = setup.D.copy() D.construct(setup.func) - BO = D.bond_order(method + option) - assert isinstance(BO, SparseAtom) + BO = D.bond_order(method + option, projection) + if projection == "atom": + assert isinstance(BO, SparseAtom) + assert BO.shape[:2] == (D.geometry.na, D.geometry.na_s) + elif projection == "orbitals": + assert isinstance(BO, SparseOrbital) + assert BO.shape[:2] == (D.geometry.no, D.geometry.no_s) def test_rho1(self, setup): D = setup.D.copy()