Total crystal angular momentum tutorial

README.md

@@ -117,7 +117,6 @@ Active contributors:
 - Miki Bonacci
 - José Castelo
 - Jorge Cervantes-Villanueva
-- Muralidhar Nalabothula
 - Riccardo Reho
 - Michele Re Fiorentin
 - Ali Esquembre-Kucukalic
@@ -131,6 +130,7 @@ Past contributors:
 - Alexandre Morlet
 - Davide Romanin
 - Daniel Murphy
+- Muralidhar Nalabothula
 
 The code is at an ongoing stage of development, help us by sending bug reports, patches and suggestions!
 

tutorial/databases_yambopy/total_crystal_angular_momentum_exciton.py

@@ -0,0 +1,123 @@
+#
+# Authors: MN
+"""
+Script to compute total crystal angular momentum
+This script loads Yambo databases, calculates the angular momentum for a specific exciton, 
+and generates a real-space visualization (.cube file).
+"""
+
+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
+
+## inputs 
+iqpt = 1
+# This is the index of the q-point (momentum transfer). 
+
+path = '.'
+# This defines the working directory where your 'GW_BSE' and 'SAVE' folders are located.
+# '.' means the current directory.
+
+gw_bse_dir = 'GW_BSE'
+# BSE Job name. This is where your bse file are store after your yambo bse run.
+
+iexe = 0
+# This is the index of the specific exciton state you want to analyze (python indexing).
+# You likely found this index by looking at the absorption spectrum or energy table previously.
+
+degen_tol = 1e-2
+# This is the energy tolerance in eV.
+# States with energy differences smaller than this value are considered "degenerate" (the same energy)
+# and will be treated together during the angular momentum mixing.
+
+# ======= Load data ========
+
+## First load the lattice db
+lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1'))
+# We load the Lattice Database (ns.db1) from the SAVE folder.
+# This file contains information about the crystal structure, unit cell geometry, and k-points.
+
+# load exciton db
+filename = 'ndb.BS_diago_Q%d' % (iqpt)
+excdb = YamboExcitonDB.from_db_file(lattice, filename=filename,
+                                    folder=os.path.join(path, gw_bse_dir), neigs = iexe + 101)
+# We load the Exciton Database.
+# Critical Note: 'neigs' (number of eigenvalues) is set to 'iexe + 100'. 
+# We load more states than the target 'iexe' to ensure we capture all states that might be degenerate 
+# with our target. If you don't load the degenerate partners, the physics will be wrong.
+
+# Load the wavefunction database
+wfdb = YamboWFDB(path=path, latdb=lattice, bands_range=[np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])])
+
+symm_mat_cart = lattice.sym_car[2] 
+# The rotation symmetry for which we want to find total_crys_angular_momentum values and simentineous eigenbasis
+# For example, here I am selecting 3th symmetry matrix which corresponds to 120 degree rotation along z (I am doing this hBN).
+# you can check this in qe scf out file when verbosity = 'high'
+# Neverthesless you can also give your custom 3x3 matrix which must satisfy Rq = q and is a symmetry of the crystal.
+
+frac_trans_cart = np.zeros(3)
+# Fractional translation vector in cart units for the corre symmetry.
+# i.e x -> Rx + t where R is symm_mat_cart and t is frac_trans_cart 
+
+sbasis = excdb.total_crys_angular_momentum(wfdb, iexe, symm_mat_cart, frac_trans_cart, degen_tol=degen_tol)
+# This is the core calculation.
+# It computes the "Total Crystal Angular Momentum" for the target exciton ('iexe').
+# It returns 'sbasis', which contains the angular momentum values and the new basis vectors.
+
+jvals = sbasis[0]
+# These are the total crystal angular momentum eigenvalues (the 'j' values) for the degenerate states.
+
+Akcv_j = sbasis[1]
+# These are the new right eigenvectors (Akcv) for the excitons i.e are simentineous eigenstates of Hbse and symmetry operator.
+# i.e they represent how to mix the original states to create states with well-defined angular momentum.
+
+print('Total crystal angular momentum : ', jvals)
+# Prints the calculated j-values to the screen.
+
+## Now we Plot read space wf
+idegen = np.array(excdb.get_degenerate(iexe + 1, eps=degen_tol), dtype=int) - 1
+# We find the indices of all states that are degenerate with our target exciton.
+# 'iexe+1' is used because some internal functions use 1-based indexing, while python uses 0-based.
+# The result 'idegen' is an array of Python indices (0-based) for these states.
+
+## Update the existing wfs with the new rotated similentious eigenbasis
+old_eig_states = excdb.Akcv[idegen].copy()
+excdb.Akcv[idegen] = Akcv_j
+# We overwrite the original exciton coefficients in the database object with our new, rotated coefficients.
+# This ensures that when we plot the state later, we are plotting the angular momentum eigenstate, 
+# not the arbitrary original state.
+
+## we break degeneracies to plot each state rather than averaging. THis is to trick the real_wf_to_cube
+## as it will always internally average the eigenstates 
+degen_tol_add = excdb.eigenvalues.max()*4
+# We prepare a large number to shift the energies.
+# Plotting tools often average degenerate states. To plot a specific 'j' state individually,
+# we need to trick the tool into thinking they have different energies.
+
+old_eig_states_ene = excdb.eigenvalues[idegen].copy()
+for i in idegen:
+    degen_tol_add += 1 
+    # Increment the shift value.
+    excdb.eigenvalues[i] += degen_tol_add
+    # We artificially shift the energy of this specific state.
+    # Now this state is energetically isolated and won't be averaged with others during plotting.
+
+# Now lets plot the real space wavefunction for each j value
+print("="*80)
+for i in range(len(jvals)):
+    print('*********** Plotting real space wf for state %d with j = %.2f ***************'
+          % (idegen[i] + 1, jvals[i]))
+    print("="*80)
+    excdb.real_wf_to_cube(iexe=idegen[i], wfdb=wfdb, fixed_postion=[0.0, 0.0, 0.0],
+                          supercell=[9,9,9], degen_tol=1e-14, wfcCutoffRy=-1,
+                          fix_particle='e', phase=True)
+    print("="*80)
+
+## restore the db back
+excdb.eigenvalues[idegen] = old_eig_states_ene
+excdb.Akcv[idegen] = old_eig_states
+
+## .cube will be dumped and use vesta to visualize it !

yambocommandline/commands/exc_sort.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 #

yambocommandline/commands/plot_exc_wf_real_space.py

@@ -1,6 +1,4 @@
 #!/usr/bin/env python3
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 import argparse

yambocommandline/commands/run_exc_irreps.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 #

yambopy/bse/exciton_irreps.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 #

yambopy/bse/exciton_matrix_elements.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 ###

yambopy/bse/exciton_spin.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Author: MN
 #

yambopy/bse/realSpace_excitonwf.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 ##

yambopy/bse/rotate_excitonwf.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 import numpy as np

yambopy/dbs/excitondb.py

@@ -26,6 +26,7 @@
 from yambopy.dbs.qpdb import YamboQPDB
 from yambopy.io.cubetools import write_cube
 from yambopy.bse.realSpace_excitonwf import ex_wf2Real
+from yambopy.bse.rotate_excitonwf import rotate_exc_wf
 
 class ExcitonList():
     """
@@ -371,6 +372,127 @@ def real_wf_to_cube(self, iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1],
                    real_wfc, sc_latvecs, atom_pos, atom_nums,
                    header='Real space exciton wavefunction')
 
+    def total_crys_angular_momentum(self, wfdb, iexe, symm_mat_cat, frac_vec_cart, degen_tol=1e-3, Dmats=None):
+        """
+        Computes the total crystal angular momentum along a given rotational symmetry axis.
+
+        This method calculates the simultaneous eigenbasis of the Hamiltonian and the
+        Symmetry operator. It verifies if the symmetry belongs to the little group of q
+        before proceeding with the representation calculations.
+
+        Parameters
+        ----------
+        wfdb : object
+            The wavefunction database object containing lattice and basis information.
+        iexe : int
+            Index of the exciton (0-based Python indexing).
+        symm_mat_cat : array_like of shape (3, 3)
+            Rotational symmetry matrix in Cartesian coordinates.
+        frac_vec_cart : array_like of shape (3,)
+            Fractional translational vector associated with the symmetry operation.
+        degen_tol : float, optional
+            Tolerance for identifying degenerate excitons in eV. Default is 1e-3.
+        Dmats : array_like, optional
+            Representation matrices for the given symmetry. If None, they are computed
+            directly from the ``wfdb``. Default is None.
+
+        Returns
+        -------
+        list
+            A list containing:
+            - w : ndarray
+                The real part of the angular momentum values (eigenvalues).
+            - sbasis_r : ndarray
+                The right symmetrized basis vectors.
+            - sbasis_l : ndarray, optional
+                The left symmetrized basis vectors. This is included only if the
+                exciton wavefunction has distinct left/right components (non-Hermitian).
+
+        Returns None
+        ------------
+        None
+            Returned if the rotation is improper (determinant < 0) or if the
+            symmetry does not belong to the little group of Q.
+        """
+        # Identify degenerate excitons
+        iexe = np.array(self.get_degenerate(iexe + 1, eps=degen_tol), dtype=int) - 1
+        print(f"Number of degenerate excitons found : {len(iexe)}. List : {iexe}")
+
+        # Check for improper rotation
+        det_r = np.linalg.det(symm_mat_cat)
+        assert det_r > 0, "Warning : Improper rotation, Only rotation are allowed."
+
+        # Prepare lattice vectors and transformation matrices
+        lat_vec = wfdb.ydb.lat
+        lat_vec_inv = np.linalg.inv(lat_vec)
+        symm_mat_red = lat_vec @ symm_mat_cat @ lat_vec_inv
+
+        # Convert the q-point to crystal coordinates
+        exc_qpt_car = self.car_qpoint
+        print(exc_qpt_car)
+        exc_qpt_crys = lat_vec @ exc_qpt_car
+        print(f"Qpt : ({exc_qpt_crys[0]:.6f}, {exc_qpt_crys[1]:.6f}, {exc_qpt_crys[2]:.6f})")
+
+        # Verify if symmetry belongs to the little group of Q
+        sq_minus_q = np.einsum('ij,j->i', symm_mat_red, exc_qpt_crys) - exc_qpt_crys
+        sq_minus_q = sq_minus_q - np.rint(sq_minus_q)
+
+        assert np.linalg.norm(sq_minus_q) < 1e-4, "Warning : The given symmetry does not belong to little group of Q."
+
+        # Compute representation matrices if not provided
+        if Dmats is None:
+            Dmats = wfdb.Dmat(
+                symm_mat=symm_mat_cat.reshape(1, 3, 3),
+                frac_vec=frac_vec_cart.reshape(1, 3),
+                time_rev=False
+            )[0]
+
+        # Retrieve exciton wavefunctions
+        # Akcv represents the coefficients of the exciton wavefunction
+        akcv = self.get_Akcv()
+        akcv_left = akcv
+
+        # Handle non-Hermitian cases where left eigenvectors differ
+        if akcv.shape[1] == 2:
+            print("Computing left ev ...")
+            akcv_left = np.linalg.inv(akcv.reshape(len(akcv), -1)).conj().T.reshape(akcv.shape)
+
+        # Select specific degenerate components
+        ak_r = akcv[iexe]
+        ak_l = akcv_left[iexe].conj()
+
+        # Calculate phase factor
+        tau_dot_k = np.exp(1j * 2 * np.pi * np.dot(self.car_qpoint, frac_vec_cart))
+
+        # Rotate the right wavefunction
+        rot_akcv = rotate_exc_wf(
+            ak_r, symm_mat_red, wfdb.kBZ, exc_qpt_crys,
+            Dmats, False, ktree=wfdb.ktree
+        )
+
+        # Compute the representation matrix
+        rep = tau_dot_k * np.einsum('m...,n...->mn', ak_l, rot_akcv, optimize=True)
+        w, v = np.linalg.eig(rep)
+
+        # Calculate rotation angle and axis via SVD
+        _, _, vt = np.linalg.svd(symm_mat_cat - np.eye(3))
+        # axis = vt[-1] / np.linalg.norm(vt[-1]) # axis is calculated but not currently used in return
+        angle = np.arccos(np.clip((np.trace(symm_mat_cat) - 1) / 2, -1.0, 1.0))
+
+        # Convert eigenvalues to angular momentum
+        w = 1j * np.log(w) / angle
+        w = w.real
+
+        # Transform basis vectors
+        sbasis_r = np.einsum('ij,i...->j...', v, ak_r, optimize=True)
+
+        if akcv.shape[1] == 2:
+            sbasis_l = np.einsum('ij,i...->j...', v, ak_l.conj(), optimize=True)
+            return [w, sbasis_r, sbasis_l]
+        else:
+            return [w, sbasis_r]
+
+
     def get_nondegenerate(self,eps=1e-4):
         """
         get a list of non-degenerate excitons
@@ -410,15 +532,15 @@ def get_sorted(self):
 
         return sort_e, sort_i 
 
-    def get_degenerate(self,index,eps=1e-4):
+    def get_degenerate(self,index,eps=1e-4,rtol=1e-3):
         """
         Get degenerate excitons
 
         Args:
             eps: maximum energy difference to consider the two excitons degenerate in eV
         """
         energy = self.eigenvalues[index-1].real
-        excitons = np.where(np.isclose(self.eigenvalues.real, energy, atol=eps))[0] + 1
+        excitons = np.where(np.isclose(self.eigenvalues.real, energy, atol=eps, rtol=rtol))[0] + 1
         return excitons.tolist()
 
     def exciton_bs(self,energies,path,excitons=(0,),debug=False):

yambopy/dbs/wfdb.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Author: MN
 

yambopy/io/cubetools.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 #

yambopy/symmetries/crystal_symmetries.py

@@ -1,5 +1,3 @@
-# Copyright (c) 2025, University of Luxembourg 
-# All rights reserved.
 #
 # Authors: MN
 #