Some progress in Excitons

pyproject.toml

@@ -13,7 +13,7 @@ authors = [
 ]
 description = "Yambopy: a pre/post-processing tool for Yambo"
 readme = "README.md"
-requires-python = ">=3.10"
+requires-python = ">=3.8"
 classifiers = [
     "Programming Language :: Python :: 3",
     "License :: OSI Approved :: GNU General Public License v2 or later (GPLv2+)",

tutorial/databases_yambopy/exc_real_wf.py

@@ -0,0 +1,40 @@
+import numpy as np
+from yambopy.dbs.excitondb import YamboExcitonDB
+from yambopy.dbs import excitondb
+from yambopy.dbs.latticedb import YamboLatticeDB
+from yambopy.dbs.wfdb import YamboWFDB
+import os
+
+iqpt = 1 # qpt index of exciton
+calc_path = '.'
+BSE_dir = 'GW_BSE'
+
+# load lattice db
+lattice = YamboLatticeDB.from_db_file(os.path.join(calc_path, 'SAVE','ns.db1'))
+
+# load exciton db
+# note in case too many excitons, load only first few with `neigs' flag
+# DO NOT forget to include all degenerate states when giving neigs flag !
+#
+filename = 'ndb.BS_diago_Q%d' % (iqpt)
+excdb = YamboExcitonDB.from_db_file(lattice, filename=filename,
+                                    folder=os.path.join(calc_path, BSE_dir),
+                                    neigs = -1)
+
+#Load the wavefunction database
+wfdb = YamboWFDB(path='.', latdb=lattice,
+                      bands_range=[np.min(excdb.table[:, 1]) - 1,
+                      np.max(excdb.table[:, 2])])
+
+## plot the exciton wavefunction with hole fixed at [0,0,0]
+# in a [1,1,1] supercell with 80 Ry wf cutoff. (give -1 to use full cutoff)
+# I want to set the degeneracy threshold to 0.01 eV
+# For example I want to plot the 3rd exciton, so iexe = 2 (python indexing )
+#
+excdb.real_wf_to_cube(iexe=2, wfdb=wfdb, fixed_postion=[0.0      ,  0.0 ,       0.0],
+                    supercell=[1,1,1], degen_tol=0.01, wfcCutoffRy=-1, fix_particle='h')
+# fixed_postion is in reduced units
+# in case, you want to plot hole density by fixing electron, set fix_particle = 'e'
+
+## .cube will be dumped and use vesta to visualize it !
+

yambopy/bse/exciton_matrix_elements.py

@@ -0,0 +1,113 @@
+###
+###
+# This file contains a genenal functions to compute
+# < S | O | S'>, where O is an operator.
+# for Ex: is O is dV_scf, then these ,matrix elements are 
+# ex-ph matrix elements. if O is S_z (spin operator), we get
+# spin matrix elements of excitons
+###
+import numpy as np
+from yambopy.kpoints import build_ktree, find_kpt
+from yambopy.tools.function_profiler import func_profile
+
+@func_profile
+def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', diagonal_only=False, ktree=None):
+    """
+    Compute the exciton matrix elements in the Tamm-Dancoff approximation.
+
+    This function calculates the matrix elements <S (final), k+q | O(q) | S' (initial), k>,
+    discarding the third term (disconnected diagram). The calculation is performed in the
+    Tamm-Dancoff approximation.
+
+    Parameters
+    ----------
+    exe_kvec : array_like
+        Exciton k-vector in crystal coordinates (k).
+    O_qvec : array_like
+        Momentum transfer vector q in crystal coordinates (q).
+    Akq : array_like
+        Wavefunction coefficients for k+q (bra wfc) with shape (n_exe_states, 1, ns, nk, nc, nv).
+    Ak : array_like
+        Wavefunction coefficients for k (ket wfc) with shape (n_exe_states, 1, ns, nk, nc, nv).
+    Omn : array_like
+        Matrix elements of the operator O in the basis of electronic states with shape (nlambda, nk, nspin, m_bnd, n_bnd).
+        ie Omn = < k+q, m, s | O(q) | n, k, s>, where m_bnd and n_bnd are final and initial bands, respectively.
+        s is spin index
+    kpts : array_like
+        K-points used to construct the BSE with shape (nk, 3).
+    contribution : str, optional
+        Specifies the contribution to include in the calculation:
+        - 'e' : Only electronic contribution.
+        - 'h' : Only hole contribution.
+        - 'b' : Both electron and hole contributions (default).
+    diagonal_only : bool, optional
+        If True, only the diagonal terms are computed. Default is False.
+
+    ktree : KDtree, optional
+        If None, will build internally, else use the user provided
+    Returns
+    -------
+    ex_O_mat : ndarray
+        The computed exciton matrix elements with shape (nlambda, n_exe_states) if diagonal_only is True,
+        or (nlambda, n_exe_states (final), n_exe_states (initial)) if diagonal_only is False.
+    """
+    # Number of arbitrary parameters (lambda) in the Omn matrix
+    nlambda = Omn.shape[0]
+    #
+    assert Akq.shape[1] == 1, "Works only with TDA."
+    # Shape of the wavefunction coefficients
+    n_exe_states, bse_calc, ns, nk, nc, nv = Akq.shape
+    #
+    # Ensure that the shapes of Akq and Ak match
+    assert Akq.shape == Ak.shape, "Wavefunction coefficient mismatch"
+    #
+    # Ensure that the contribution parameter is valid
+    assert contribution in ['b', 'e', 'h'], "Allowed values for contribution are 'b', 'e', 'h'"
+    #
+    # Build a k-point tree for efficient k-point searching
+    if ktree is None : ktree = build_ktree(kpts)
+    #
+    # Find the indices of k-Q-q and k-q in the k-point tree
+    idx_k_minus_Q_minus_q = find_kpt(ktree, kpts - O_qvec[None, :] - exe_kvec[None, :])  # k-Q-q
+    idx_k_minus_q = find_kpt(ktree, kpts - O_qvec[None, :])  # k-q
+    #
+    # Extract the occupied and unoccupied parts of the Omn matrix
+    Occ = Omn[:, idx_k_minus_q, :, nv:, nv:].transpose(0,2,1,3,4)  # Occupied part
+    Ovv = Omn[:, idx_k_minus_Q_minus_q, :, :nv, :nv].transpose(0,2,1,3,4)  # conduction part
+    #
+    # Ensure the arrays are C-contiguous to reduce cache misses
+    Ak_electron = np.ascontiguousarray(Ak[:,0][:,:,idx_k_minus_q, ...])
+    Akq_conj = Akq[:,0].reshape(n_exe_states, -1).conj()
+    #
+    # Initialize the output matrix
+    if diagonal_only:
+        ex_O_mat = np.zeros((nlambda, n_exe_states), dtype=Ak.dtype)  # (nlambda, final, initial)
+    else:
+        ex_O_mat = np.zeros((nlambda, n_exe_states, n_exe_states), dtype=Ak.dtype)  # (nlambda, final, initial)
+    #
+    # Loop over the arbitrary parameters (lambda)
+    for il in range(nlambda):
+        # Compute the electron contribution
+        if contribution == 'e' or contribution == 'b':
+            tmp_wfc = Occ[il][None, ...] @ Ak_electron
+        #
+        # Compute the hole contribution and subtract from the electron contribution
+        if contribution == 'h' or contribution == 'b':
+            tmp_h = -Ak[:,0] @ Ovv[il][None, ...]
+            if contribution == 'b':
+                tmp_wfc += tmp_h
+            else:
+                tmp_wfc = tmp_h
+        #
+        # Reshape the temporary wavefunction coefficients
+        tmp_wfc = tmp_wfc.reshape(n_exe_states, -1)
+        #
+        # Compute the final matrix elements
+        if diagonal_only:
+            ex_O_mat[il] = np.sum(Akq_conj * tmp_wfc, axis=-1)
+        else:
+            np.matmul(Akq_conj, tmp_wfc.T, out=ex_O_mat[il])
+    #
+    # Return the computed exciton matrix elements
+    return ex_O_mat
+

yambopy/bse/exciton_spin.py

@@ -0,0 +1,179 @@
+import os
+import numpy as np
+from yambopy.dbs.excitondb import YamboExcitonDB
+from yambopy.dbs.latticedb import YamboLatticeDB
+from yambopy.dbs.wfdb import YamboWFDB
+from .exciton_matrix_elements import exciton_X_matelem
+from yambopy.tools.degeneracy_finder import find_degeneracy_evs
+
+
+def compute_exciton_spin(lattice, excdb, wfdb, elec_sz, contribution='b',diagonal=False):
+    """
+    Compute the spin matrix elements <S'|S_z|S> for excitons.
+
+    This function calculates the spin matrix elements for excitons using the
+    wavefunctions and spin operators. The spin matrix is computed in the basis
+    of exciton states, and off-diagonal elements are included. Diagonalization
+    of the matrix in degenerate subspaces is required to obtain the spin values.
+
+    Parameters
+    ----------
+    lattice : latticedb
+        Lattice database
+    excdb : exciton db
+        Exciton Database 
+    wfdb : wfc db
+        wavefunction Database 
+    elec_sz : ndarray
+        Electron spin matrix elements (nk, nbnds. nbnds).
+    contribution : str, optional
+        Specifies which contribution to compute:
+        - 'b': Total spin (default).
+        - 'e': Electron spin only.
+        - 'h': Hole spin only.
+    diagonal : bool, optional
+        If True, only diagonal spin elements are computed. Default is False.
+
+    Returns
+    -------
+    exe_Sz : ndarray
+        Spin matrix elements for excitons with shape (nstates, nstates).
+    """
+    #
+    # Ensure the calculation is valid only for spinor wavefunctions
+    assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2"
+    #
+    # Sanity check
+    assert np.min(excdb.table[:, 1]) - 1 == wfdb.min_bnd, \
+            "wfdb and exciton db are inconsistant (Bands)"
+    ## sanity check
+    assert np.max(excdb.table[:, 2]) == wfdb.min_bnd + wfdb.nbands, \
+            "wfdb and exciton db are inconsistant (Bands)"
+    #
+    assert elec_sz.shape == (wfdb.nkBZ, wfdb.nbands, wfdb.nbands)
+    # get Akcv
+    Akcv = excdb.get_Akcv()
+    #
+    # Get the exciton q-point in Cartesian coordinates
+    excQpt = excdb.car_qpoint
+    #
+    # Convert the q-point to crystal coordinates
+    excQpt = lattice.lat @ excQpt
+    #
+    # Compute the exciton spin matrix elements <S'|S_z|S>
+    exe_Sz = exciton_X_matelem(excQpt, np.array([0, 0, 0]), Akcv,
+                               Akcv, elec_sz[None,:,None,...], wfdb.kBZ,
+                               diagonal_only=diagonal,contribution=contribution)
+    #
+    return exe_Sz[0]
+
+
+
+
+
+def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1,
+                          nstates=-1, contribution='b', degen_tol = 1e-2,
+                          sz=0.5 * np.array([[1, 0], [0, -1]]),
+                          return_dbs_and_spin=True):
+    """
+
+
+    Description
+    -----------
+    Compute expectation value of S_z operator for excitons.
+
+    Parameters
+    ----------
+    path : str, optional
+        Path to the directory containing calculation SAVE and BSE folder.
+        Default: '.' (current directory)
+    bse_dir : str, optional
+        Directory containing BSE calculation data. Default: 'SAVE'
+    iqpt : int or array-like, optional
+        Q-point index or list of Q-point indices to analyze. Default: 1
+        (Fortran indexing)
+    nstates : int, optional
+        Number of excitonic states to consider. Use -1 for all states. Default: -1
+    contribution : str, optional
+        Which contribution to compute:
+        - 'b': both electron and hole (default)
+        - 'e': electron only
+        - 'h': hole only
+    degen_tol : float, optional
+        Tolerance for detecting degenerate states. Default: 1e-2
+    sz : ndarray, optional
+        S_z operator matrix representation. Default: 0.5 * np.array([[1, 0], [0, -1]])
+    return_dbs_and_spin : bool, optional
+        If True, returns both spin values and database objects. Default: True
+
+    Returns
+    -------
+    exe_Sz : ndarray
+        Array containing S_z expectation values for excitonic states
+    dbs_objects : list, optional
+        If return_dbs_and_spin=True, returns [lattice, wfdb, excdb, elec_sz] database objects
+    Examples
+    --------
+    Compute the total spin matrix elements for excitons:
+
+    >>> import numpy as np
+    >>> from yambopy.bse.exciton_spin import compute_exc_spin_iqpt
+    >>> Sz_exe = compute_exc_spin_iqpt(bse_dir='GW_BSE', nstates=2)
+    >>> print(Sz_exe)
+
+    Compute only the electron spin contribution:
+
+    >>> Sz_exe = compute_exc_spin_iqpt(bse_dir='GW_BSE', nstates=2, contribution='e')
+
+    Compute only the hole spin contribution:
+
+    >>> Sz_exe = compute_exc_spin_iqpt(bse_dir='GW_BSE', nstates=2, contribution='h')
+    """
+    #
+    ## Check if it single Q or multiple Q's
+    if np.isscalar(iqpt): iqpt = [iqpt]
+    else : iqpt = list(iqpt)
+    # Load the lattice database
+    lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1'))
+    ## load exbds
+    excdb = []
+    for iq in iqpt:
+        filename = 'ndb.BS_diago_Q%d' % (iq)
+        excdb.append(YamboExcitonDB.from_db_file(lattice, filename=filename,
+                                                 folder=os.path.join(path, bse_dir),
+                                                 Load_WF=True, neigs=nstates))
+    # Load the wavefunction database
+    wfdb = YamboWFDB(path=path, latdb=lattice,
+                      bands_range=[np.min(excdb[0].table[:, 1]) - 1,
+                    np.max(excdb[0].table[:, 2])])
+    #
+    # Compute the spin matrix elements in the BZ
+    elec_sz = wfdb.get_spin_m_e_BZ(s_z=sz)
+    #
+    exe_Sz = []
+    for ixdb in excdb:
+        smat = compute_exciton_spin(lattice, ixdb,
+                                           wfdb, elec_sz,
+                                           contribution=contribution,
+                                           diagonal=False)
+        smat = get_spinvals(smat, ixdb.eigenvalues, atol=degen_tol)
+        ss_tmp = []
+        for i in smat: ss_tmp = ss_tmp + list(i)
+        exe_Sz.append(ss_tmp)
+    #
+    exe_Sz = np.array(exe_Sz)
+    if return_dbs_and_spin : return exe_Sz,[lattice, wfdb, excdb, elec_sz]
+    else : return exe_Sz
+
+
+
+def get_spinvals(spin_matrix, eigenvalues, atol=1e-3, rtol=1e-3):
+    degen_idx = find_degeneracy_evs(eigenvalues,atol=atol, rtol=rtol)
+    spins = []
+    for id in degen_idx:
+        w = np.linalg.eigvals(spin_matrix[id,:][:,id])
+        spins.append(w)
+    return spins
+
+
+