From 3fe946ca2e958d8c5fc6ff1774b85c4c3eab4292 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 16:28:22 +0100 Subject: [PATCH 01/73] Add exciton phonon routines --- yambopy/bse/exciton_matrix_elements.py | 108 +++++++++++++++++++++++++ yambopy/bse/rotate_excitonwf.py | 62 ++++++++++++++ 2 files changed, 170 insertions(+) create mode 100644 yambopy/bse/exciton_matrix_elements.py create mode 100644 yambopy/bse/rotate_excitonwf.py diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py new file mode 100644 index 00000000..fb0be3b8 --- /dev/null +++ b/yambopy/bse/exciton_matrix_elements.py @@ -0,0 +1,108 @@ +### +### +# 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 + +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 , + 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. + O_qvec : array_like + Momentum transfer vector q in crystal coordinates. + Akq : array_like + Wavefunction coefficients for k+q (bra wfc) with shape (n_exe_states, nk, nc, nv). + Ak : array_like + Wavefunction coefficients for k (ket wfc) with shape (n_exe_states, nk, nc, nv). + Omn : array_like + Matrix elements of the operator O in the basis of electronic states with shape (nlambda, nk, nm, nn). + 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] + # + # Shape of the wavefunction coefficients + n_exe_states, 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 not ktree : 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:] # Occupied part + Ovv = Omn[:, idx_k_minus_Q_minus_q, :nv, :nv] # conduction part + # + # Ensure the arrays are C-contiguous to reduce cache misses + Ak_electron = np.ascontiguousarray(Ak[:, idx_k_minus_q, ...]) + Akq_conj = Akq.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 @ 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 + diff --git a/yambopy/bse/rotate_excitonwf.py b/yambopy/bse/rotate_excitonwf.py new file mode 100644 index 00000000..54b8b620 --- /dev/null +++ b/yambopy/bse/rotate_excitonwf.py @@ -0,0 +1,62 @@ +import numpy as np +from yambopy.kpoints import build_ktree, find_kpt + +def rotate_exc_Ak(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=None): + """ + Rotate the exciton wavefunction Ak using symmetry operations. + + This function applies a symmetry operation to the exciton wavefunction Ak, + which is represented in the basis of electronic states. The rotation is + performed using the symmetry matrix in reduced coordinates and the + corresponding representation matrices. + + Parameters + ---------- + Ak : array_like + Exciton wavefunction coefficients with shape (n_exe_states, nk, nc, nv). + symm_mat_red : array_like + Symmetry matrix in reduced coordinates with shape (3, 3). + kpoints : array_like + K-points in the full Brillouin zone (crystal coordinates) with shape (nk, 3). + exe_qpt : array_like + Momentum of the exciton (q-point) in crystal coordinates with shape (3,). + dmats : array_like + Representation matrices for the symmetry operation with shape (nk, Rk_band, k_band). + time_rev : bool + If True, apply time-reversal symmetry to the wavefunction. + ktree : object, optional + Pre-built k-point tree for efficient k-point searching. If not provided, it will be built. + + Returns + ------- + rot_Ak : ndarray + Rotated exciton wavefunction coefficients with the same shape as Ak. + """ + # Initialize the rotated Ak array + rot_Ak = np.zeros(Ak.shape, dtype=Ak.dtype) + nk, nc, nv = Ak.shape[1:] + + # Build a k-point tree if not provided + if not ktree: + ktree = build_ktree(kpoints) + + # Compute the indices of Rk and Rk - q + Rkpts = kpoints @ symm_mat_red.T # Rotated k-points + k_minus_q = kpoints - exe_qpt[None, :] # k - q + idx_Rk = find_kpt(ktree, Rkpts) # Indices of rotated k-points + idx_k_minus_q = find_kpt(ktree, k_minus_q) # Indices of k - q + + # Extract the conduction and valence parts of the representation matrices + Dcc = dmats[:, nv:, nv:] # Conduction band part + Dvv = dmats[idx_k_minus_q, :nv, :nv].conj() # Valence band part (conjugated) + + # Apply time-reversal symmetry if required + Ak_tmp = Ak + if time_rev: + Ak_tmp = Ak.conj() + + # Rotate the Ak wavefunction using the representation matrices + rot_Ak[:, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp) @ (Dvv.transpose(0, 2, 1)[None, ...]) + + return rot_Ak + From 4b5305c137af4be4bcd00f6a5acc1e13a08546a4 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 17:36:22 +0100 Subject: [PATCH 02/73] Add option to load only subset of exe wfs --- yambopy/bse/rotate_excitonwf.py | 2 +- yambopy/dbs/excitondb.py | 66 ++++++++++++++++++++++++++++----- 2 files changed, 57 insertions(+), 11 deletions(-) diff --git a/yambopy/bse/rotate_excitonwf.py b/yambopy/bse/rotate_excitonwf.py index 54b8b620..0e9ca918 100644 --- a/yambopy/bse/rotate_excitonwf.py +++ b/yambopy/bse/rotate_excitonwf.py @@ -1,7 +1,7 @@ import numpy as np from yambopy.kpoints import build_ktree, find_kpt -def rotate_exc_Ak(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=None): +def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=None): """ Rotate the exciton wavefunction Ak using symmetry operations. diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index b8c6ad50..c45af278 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -82,13 +82,15 @@ def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,spin_pol='no',ca self.table = table self.eigenvectors = eigenvectors self.spin_pol = spin_pol + self.neigs = len(eigenvalues) @classmethod - def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True): + def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True, neigs=-1): """ Initialize this class from a file Set `Read_WF=False` to avoid reading eigenvectors for faster IO and memory efficiency. + If neigs < 0 ; all eigen values (vectors) are loaded or else first neigs are loaded """ path_filename = os.path.join(folder,filename) if not os.path.isfile(path_filename): @@ -98,6 +100,13 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True) Qpt = filename.split("Q",1)[1] with Dataset(path_filename) as database: + #energies + eig = database.variables['BS_Energies'][:]*ha2ev + eigenvalues = eig[:,0]+eig[:,1]*I + neig_full = len(eigenvalues) + if neigs < 0 or neigs > neig_full: neigs = neig_full + eigenvalues = eigenvalues[:neigs] + if 'BS_left_Residuals' in list(database.variables.keys()): # MN: using complex views instead of a+I*b copies to avoid memory duplication # Old (yet instructive) memory duplication code @@ -105,13 +114,13 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True) #rer,imr = database.variables['BS_right_Residuals'][:].T #l_residual = rel+iml*I #r_residual = rer+imr*I - l_residual = database.variables['BS_left_Residuals'][:] - r_residual = database.variables['BS_right_Residuals'][:] + l_residual = database['BS_left_Residuals'][:neigs,...].data + r_residual = database['BS_right_Residuals'][:neigs,...].data l_residual = l_residual.view(dtype=CmplxType(l_residual)).reshape(len(l_residual)) r_residual = r_residual.view(dtype=CmplxType(r_residual)).reshape(len(r_residual)) if 'BS_Residuals' in list(database.variables.keys()): # Compatibility with older Yambo versions - rel,iml,rer,imr = database.variables['BS_Residuals'][:].T + rel,iml,rer,imr = database['BS_Residuals'][:neigs,...].data.T l_residual = rel+iml*I r_residual = rer+imr*I @@ -121,19 +130,15 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True) car_qpoint = database.variables['Q-point'][:]/lattice.alat if Qpt=="1": car_qpoint = np.zeros(3) - #energies - eig = database.variables['BS_Energies'][:]*ha2ev - eigenvalues = eig[:,0]+eig[:,1]*I - #eigenvectors table = None eigenvectors = None if Load_WF and 'BS_EIGENSTATES' in database.variables: - eiv = database.variables['BS_EIGENSTATES'][:] + eiv = database['BS_EIGENSTATES'][:neigs,...].data #eiv = eiv[:,:,0] + eiv[:,:,1]*I #eigenvectors = eiv eigenvectors = eiv.view(dtype=CmplxType(eiv)).reshape(eiv.shape[:-1]) - table = database.variables['BS_TABLE'][:].T.astype(int) + table = np.rint(database.variables['BS_TABLE'][:].T).astype(int) spin_vars = [int(database.variables['SPIN_VARS'][:][0]), int(database.variables['SPIN_VARS'][:][1])] if spin_vars[0] == 2 and spin_vars[1] == 1: @@ -230,6 +235,47 @@ def write_sorted(self,prefix='yambo'): for i,n in sort_i: f.write("%3d %12.8lf %12.8e\n"%(n+1,eig[n],i)) + def convert_to_kcv(self): + """ + Convert eigenvectors from (neigs,BS_table) -> (neigs,k,c,v) + For now, only works for nspin = 1. nspinor = 1/2 also works. + """ + assert self.spin_pol == 'no', "Rearrange_Akcv works only for nspin = 1" + # + if self.eigenvectors is None: return None + eig_wfcs = self.eigenvectors + # + nk = self.nkpoints + nv = self.nvbands + nc = self.ncbands + # Make sure nc * nv * nk = BS_TABLE length + table_len = nk*nv*nc + assert table_len == self.table.shape[0], "BS_TABLE length not equal to nc * nv * nk" + + assert eig_wfcs.shape[-1]//table_len == 1, "rearranged_Akcv works only in TDA" + # + v_min = np.min(self.table[:,1]) + c_min = np.min(self.table[:,2]) + bs_table0 = self.table[:,0]-1 + bs_table1 = self.table[:,1] - v_min + bs_table2 = self.table[:,2] - c_min + # + eig_wfcs_returned = np.zeros(eig_wfcs.shape,dtype=eig_wfcs.dtype) + # + sort_idx = bs_table0*nc*nv + bs_table2*nv + bs_table1 + # + eig_wfcs_returned[:,sort_idx] = eig_wfcs[...,:table_len] + eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) + # + # check if this is coupling . + # if eig_wfcs.shape[-1]//table_len == 2: + # eig_wfcs_returned[:,sort_idx+table_len] = eig_wfcs[...,table_len:] + # eig_wfcs_returned = eig_wfcs_returned.reshape(-1,2,nk,nc,nv) + # else : + # eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) + # + self.eigenvectors = eig_wfcs_returned + def get_nondegenerate(self,eps=1e-4): """ get a list of non-degenerate excitons From 7a535c984cb2356ed86ca1bb06d61f4b94b38862 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 18:53:09 +0100 Subject: [PATCH 03/73] Exciton spin --- yambopy/bse/exciton_spin.py | 37 +++++++++++++++++++++++++++++++++++++ yambopy/dbs/excitondb.py | 1 - 2 files changed, 37 insertions(+), 1 deletion(-) create mode 100644 yambopy/bse/exciton_spin.py diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py new file mode 100644 index 00000000..648ccda5 --- /dev/null +++ b/yambopy/bse/exciton_spin.py @@ -0,0 +1,37 @@ +## Compute diagonal expectation value of spin operator +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 + +def compute_exciton_spin(path='.',bse_dir='',iqpt=1, nstates=-1, + sz = 0.5*np.array([[1, 0], [0, -1]]) ): + ## Computing + filename = 'ndb.BS_diago_Q%d'%(iqpt) + # + lattice = YamboLatticeDB.from_db_file(os.path.join(path,'SAVE','ns.db1')) + # + excdb = YamboExcitonDB.from_db_file(lattice,filename='ndb.BS_diago_Q1', + folder=os.path.join(path,bse_dir), + Load_WF=True, neigs=nstates) + # + wfdb = YamboWFDB(path=path,bands_range=[np.min(excdb.table[:,1])-1, + np.max(excdb.table[:,2])] ) + # + print(np.min(excdb.table[:,1])-1,np.max(excdb.table[:,2])) + assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" + # + excdb.convert_to_kcv() + # + elec_sz = wfdb.get_spin_m_e_BZ(self,s_z=sz) + # + excQpt = excdb.Qpt + # convert to crystal units + excQpt = lattice.lat@excQpt + # + exe_Sz = exciton_X_matelem(excQpt, np.array([0,0,0]), excdb.eigenvectors, + excdb.eigenvectors, elec_sz[None,...], wfdb.kBZ, + diagonal_only=True) + print(exe_Sz) diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index c45af278..fe9722a2 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -82,7 +82,6 @@ def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,spin_pol='no',ca self.table = table self.eigenvectors = eigenvectors self.spin_pol = spin_pol - self.neigs = len(eigenvalues) @classmethod def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True, neigs=-1): From 88d5689cae13d23f746fbb0a4bb9ab2a02ce4e52 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 19:18:27 +0100 Subject: [PATCH 04/73] fix bug in wfdb --- yambopy/bse/exciton_spin.py | 6 +++--- yambopy/dbs/wfdb.py | 6 +++++- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 648ccda5..e73f285e 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -6,7 +6,7 @@ from yambopy.dbs.wfdb import YamboWFDB from .exciton_matrix_elements import exciton_X_matelem -def compute_exciton_spin(path='.',bse_dir='',iqpt=1, nstates=-1, +def compute_exciton_spin(path='.',bse_dir='SAVE',iqpt=1, nstates=-1, sz = 0.5*np.array([[1, 0], [0, -1]]) ): ## Computing filename = 'ndb.BS_diago_Q%d'%(iqpt) @@ -25,9 +25,9 @@ def compute_exciton_spin(path='.',bse_dir='',iqpt=1, nstates=-1, # excdb.convert_to_kcv() # - elec_sz = wfdb.get_spin_m_e_BZ(self,s_z=sz) + elec_sz = wfdb.get_spin_m_e_BZ(s_z=sz) # - excQpt = excdb.Qpt + excQpt = excdb.car_qpoint # convert to crystal units excQpt = lattice.lat@excQpt # diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 72951452..7206a2d6 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -97,9 +97,13 @@ def read(self, bands_range=[]): lat_param = self.ydb.alat # K-points in iBZ (crystal units) - self.kpts_iBZ = self.ydb.iku_kpoints / lat_param[None, :] + self.kpts_iBZ = self.ydb.ibz_kpoints / lat_param[None, :] self.kpts_iBZ = self.kpts_iBZ @ lat_vec + # K-points in BZ (crystal units) + self.kBZ = self.ydb.iku_kpoints / lat_param[None, :] + self.kBZ = self.kpts_iBZ @ lat_vec + # G-vectors in cartesian units G_vec = ns_db1['G-VECTORS'][...].data.T wfc_grid = ns_db1['WFC_GRID'][...].data From 9b0dee093f8f19521a0b79ef0feb10eac15b067b Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 19:29:45 +0100 Subject: [PATCH 05/73] fix bug in wfdb --- yambopy/dbs/wfdb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 7206a2d6..e0c499b5 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -102,7 +102,7 @@ def read(self, bands_range=[]): # K-points in BZ (crystal units) self.kBZ = self.ydb.iku_kpoints / lat_param[None, :] - self.kBZ = self.kpts_iBZ @ lat_vec + self.kBZ = self.kBZ @ lat_vec # G-vectors in cartesian units G_vec = ns_db1['G-VECTORS'][...].data.T From 20d698567a49e3794e4b8565973f7dafeebe2767 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 19:55:31 +0100 Subject: [PATCH 06/73] exciton spin working --- yambopy/bse/exciton_spin.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index e73f285e..f2a65b94 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -1,4 +1,3 @@ -## Compute diagonal expectation value of spin operator import os import numpy as np from yambopy.dbs.excitondb import YamboExcitonDB @@ -8,7 +7,9 @@ def compute_exciton_spin(path='.',bse_dir='SAVE',iqpt=1, nstates=-1, sz = 0.5*np.array([[1, 0], [0, -1]]) ): - ## Computing + ## Computing + # Note that it also computes off-diagonal. + ## One needs to diagonalize in degnerate subspace to get the spin of excitons filename = 'ndb.BS_diago_Q%d'%(iqpt) # lattice = YamboLatticeDB.from_db_file(os.path.join(path,'SAVE','ns.db1')) @@ -33,5 +34,6 @@ def compute_exciton_spin(path='.',bse_dir='SAVE',iqpt=1, nstates=-1, # exe_Sz = exciton_X_matelem(excQpt, np.array([0,0,0]), excdb.eigenvectors, excdb.eigenvectors, elec_sz[None,...], wfdb.kBZ, - diagonal_only=True) - print(exe_Sz) + diagonal_only=False) + print("Note : This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") + return exe_Sz From 52bdaa99b225ce599465d54fc7225b12e8e72c4f Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 20:19:31 +0100 Subject: [PATCH 07/73] Exciton spin example --- yambopy/bse/exciton_spin.py | 97 ++++++++++++++++++++++++++++--------- 1 file changed, 73 insertions(+), 24 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index f2a65b94..1c640d12 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -5,35 +5,84 @@ from yambopy.dbs.wfdb import YamboWFDB from .exciton_matrix_elements import exciton_X_matelem -def compute_exciton_spin(path='.',bse_dir='SAVE',iqpt=1, nstates=-1, - sz = 0.5*np.array([[1, 0], [0, -1]]) ): - ## Computing - # Note that it also computes off-diagonal. - ## One needs to diagonalize in degnerate subspace to get the spin of excitons - filename = 'ndb.BS_diago_Q%d'%(iqpt) - # - lattice = YamboLatticeDB.from_db_file(os.path.join(path,'SAVE','ns.db1')) +def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, + sz=0.5 * np.array([[1, 0], [0, -1]])): + """ + Compute the spin matrix elements 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 + ---------- + path : str, optional + Path to the directory containing the SAVE folder. Default is '.'. + bse_dir : str, optional + Directory containing the BSE data. Default is 'SAVE'. + iqpt : int, optional + Index of the q-point for which the exciton spin is computed. Default is 1. + nstates : int, optional + Number of exciton states to consider. If -1, all states are included. Default is -1. + sz : array_like, optional + Spin-z operator matrix in the basis of spinor wavefunctions. Default is 0.5 * [[1, 0], [0, -1]]. + + Returns + ------- + exe_Sz : ndarray + Spin matrix elements for excitons with shape (nstates, nstates). + + -------------------------------------------------------------------- + Example: + -------------------------------------------------------------------- + import numpy as np + from yambopy.bse.exciton_spin import compute_exciton_spin # - excdb = YamboExcitonDB.from_db_file(lattice,filename='ndb.BS_diago_Q1', - folder=os.path.join(path,bse_dir), + Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4) + w = np.linalg.eigvals(Sz_exe) + print(w) ## spin values of excitons + """ + # Filename for the BSE diagonalization data + filename = 'ndb.BS_diago_Q%d' % (iqpt) + + # Load the lattice database + lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) + + # Load the exciton database + excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, + folder=os.path.join(path, bse_dir), Load_WF=True, neigs=nstates) - # - wfdb = YamboWFDB(path=path,bands_range=[np.min(excdb.table[:,1])-1, - np.max(excdb.table[:,2])] ) - # - print(np.min(excdb.table[:,1])-1,np.max(excdb.table[:,2])) + + # Load the wavefunction database + wfdb = YamboWFDB(path=path, bands_range=[np.min(excdb.table[:, 1]) - 1, + np.max(excdb.table[:, 2])]) + + # Print the band range for debugging + print(np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])) + + # Ensure the calculation is valid only for spinor wavefunctions assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" - # + + # Convert the exciton database to kcv format excdb.convert_to_kcv() - # + + # Compute the spin matrix elements in the BZ elec_sz = wfdb.get_spin_m_e_BZ(s_z=sz) - # + + # Get the exciton q-point in Cartesian coordinates excQpt = excdb.car_qpoint - # convert to crystal units - excQpt = lattice.lat@excQpt - # - exe_Sz = exciton_X_matelem(excQpt, np.array([0,0,0]), excdb.eigenvectors, - excdb.eigenvectors, elec_sz[None,...], wfdb.kBZ, + + # Convert the q-point to crystal coordinates + excQpt = lattice.lat @ excQpt + + # Compute the exciton spin matrix elements + exe_Sz = exciton_X_matelem(excQpt, np.array([0, 0, 0]), excdb.eigenvectors, + excdb.eigenvectors, elec_sz[None, ...], wfdb.kBZ, diagonal_only=False) - print("Note : This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") + + # Print a note about the spin matrix + print("Note: This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") + return exe_Sz + From 728a197a8c1d78937c82a17991a45552a7acd1e1 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 21:28:37 +0100 Subject: [PATCH 08/73] remove debug comment in exciton spin --- yambopy/bse/exciton_spin.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 1c640d12..4d9eab58 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -58,9 +58,6 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, wfdb = YamboWFDB(path=path, bands_range=[np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])]) - # Print the band range for debugging - print(np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])) - # Ensure the calculation is valid only for spinor wavefunctions assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" From 62c6a9b93105a9363f9881816d4a1f6f3ab20ecf Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 22:14:56 +0100 Subject: [PATCH 09/73] Add option to compute contribution --- yambopy/bse/exciton_spin.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 4d9eab58..48930aa5 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -5,7 +5,7 @@ from yambopy.dbs.wfdb import YamboWFDB from .exciton_matrix_elements import exciton_X_matelem -def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, +def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribution='b', sz=0.5 * np.array([[1, 0], [0, -1]])): """ Compute the spin matrix elements for excitons. @@ -27,6 +27,9 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, Number of exciton states to consider. If -1, all states are included. Default is -1. sz : array_like, optional Spin-z operator matrix in the basis of spinor wavefunctions. Default is 0.5 * [[1, 0], [0, -1]]. + contribution : char, optional + 'b','e', 'h'. If 'b' total spin is computed. 'e'/'h' for only electron/hole spin. + Default is 'b' (total spin) Returns ------- @@ -39,9 +42,20 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, import numpy as np from yambopy.bse.exciton_spin import compute_exciton_spin # + # + # Total spin Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4) + # + # Only Electron spin + Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4,contribution='e') + # + # Only Hole spin + Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4,contribution='h') + # w = np.linalg.eigvals(Sz_exe) print(w) ## spin values of excitons + # + # """ # Filename for the BSE diagonalization data filename = 'ndb.BS_diago_Q%d' % (iqpt) @@ -73,13 +87,13 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, # Convert the q-point to crystal coordinates excQpt = lattice.lat @ excQpt - # Compute the exciton spin matrix elements + # Compute the exciton spin matrix elements exe_Sz = exciton_X_matelem(excQpt, np.array([0, 0, 0]), excdb.eigenvectors, excdb.eigenvectors, elec_sz[None, ...], wfdb.kBZ, - diagonal_only=False) + diagonal_only=False,contribution=contribution) # Print a note about the spin matrix - print("Note: This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") + #print("Note: This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") return exe_Sz From 47fba5da96e79ac934406561aeb1a14357f68e47 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 23:13:34 +0100 Subject: [PATCH 10/73] fix bugs in exciton spin --- yambopy/bse/exciton_spin.py | 31 +++++++++++++++++++++++-------- yambopy/dbs/excitondb.py | 5 +++++ yambopy/dbs/wfdb.py | 13 +++++++++---- 3 files changed, 37 insertions(+), 12 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 48930aa5..642f05a6 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -6,7 +6,7 @@ from .exciton_matrix_elements import exciton_X_matelem def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribution='b', - sz=0.5 * np.array([[1, 0], [0, -1]])): + sz=0.5 * np.array([[1, 0], [0, -1]]), lattice=None, wfdb=None, excdb=None): """ Compute the spin matrix elements for excitons. @@ -31,6 +31,8 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut 'b','e', 'h'. If 'b' total spin is computed. 'e'/'h' for only electron/hole spin. Default is 'b' (total spin) + if latticedb, wfdb, excitondb are not None, then the provided db is used and is not read again + # Returns ------- exe_Sz : ndarray @@ -44,14 +46,16 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut # # # Total spin - Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4) + Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=2) # # Only Electron spin - Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4,contribution='e') + Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=2,contribution='e') # # Only Hole spin - Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=4,contribution='h') + Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=2,contribution='h') + # # + # The first two states are degenerate so, we diagonalize in 2x2 sub block w = np.linalg.eigvals(Sz_exe) print(w) ## spin values of excitons # @@ -61,17 +65,28 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut filename = 'ndb.BS_diago_Q%d' % (iqpt) # Load the lattice database - lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) + if not lattice: + lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) # Load the exciton database - excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, + if not excdb: + excdb = 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, bands_range=[np.min(excdb.table[:, 1]) - 1, - np.max(excdb.table[:, 2])]) + if not wfdb: + wfdb = YamboWFDB(path=path, latdb=lattice, + bands_range=[np.min(excdb.table[:, 1]) - 1, + np.max(excdb.table[:, 2])]) + ## sanity check + assert np.min(excdb.table[:, 1]) - 1 == wfdb.min_bnd, \ + "wfdb and exciton db are inconsistant (Bands)" + # + assert np.max(excdb.table[:, 2]) == wfdb.min_bnd + wfdb.nbands, \ + "wfdb and exciton db are inconsistant (Bands)" + # # Ensure the calculation is valid only for spinor wavefunctions assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index fe9722a2..f00d1ded 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -82,6 +82,9 @@ def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,spin_pol='no',ca self.table = table self.eigenvectors = eigenvectors self.spin_pol = spin_pol + self.form = 'bs_tab' ## ('bs_tab or 'kcv') + # if form is bs_tab, them eigenvectors are stored in (nexe,BS_table) + # if form == 'kcv', they are stored as (nexe, nk, nc, nv) @classmethod def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True, neigs=-1): @@ -239,6 +242,8 @@ def convert_to_kcv(self): Convert eigenvectors from (neigs,BS_table) -> (neigs,k,c,v) For now, only works for nspin = 1. nspinor = 1/2 also works. """ + if self.form == 'kcv': return + self.form = 'kcv' assert self.spin_pol == 'no', "Rearrange_Akcv works only for nspin = 1" # if self.eigenvectors is None: return None diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index e0c499b5..16741e13 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -55,7 +55,7 @@ class YamboWFDB: to_real_space(wfc_tmp, gvec_tmp, grid=[]): Convert wavefunctions to real space. """ - def __init__(self, path=None, save='SAVE', filename='ns.wf', bands_range=[]): + def __init__(self, path=None, save='SAVE', filename='ns.wf', bands_range=[], latdb=None): """ Initialize the YamboWFDB class. @@ -72,14 +72,15 @@ def __init__(self, path=None, save='SAVE', filename='ns.wf', bands_range=[]): self.wf_expanded = False ## if true, then wfc's are expanded over full BZ # Read wavefunctions - self.read(bands_range=bands_range) + self.read(bands_range=bands_range, latdb=latdb) - def read(self, bands_range=[]): + def read(self, bands_range=[], latdb=None): """ Read wavefunctions from the file. Args: bands_range (list, optional): Range of bands to load. Defaults to all bands. + latdb : latticedb, if None (default), it will be created internally """ path = self.path filename = self.filename @@ -87,7 +88,11 @@ def read(self, bands_range=[]): # Open the ns.db1 database to get essential data try: ns_db1_fname = os.path.join(path, 'ns.db1') - self.ydb = YamboLatticeDB.from_db_file(ns_db1_fname, Expand=True) + if latdb : + if not hasattr(latdb,'ibz_kpoints'): latdb.expand_kpoints() + self.ydb = latdb + else : + self.ydb = YamboLatticeDB.from_db_file(ns_db1_fname, Expand=True) ## total kpoints in full BZ self.nkBZ = len(self.ydb.symmetry_indexes) # From 634f579b5b20e360e0d55f8004d8d91e18cedb81 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 23:51:01 +0100 Subject: [PATCH 11/73] Support multiple q's given at a time exciton spin --- yambopy/bse/exciton_spin.py | 77 +++++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 33 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 642f05a6..a26e7ecf 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -6,7 +6,8 @@ from .exciton_matrix_elements import exciton_X_matelem def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribution='b', - sz=0.5 * np.array([[1, 0], [0, -1]]), lattice=None, wfdb=None, excdb=None): + sz=0.5 * np.array([[1, 0], [0, -1]]), + lattice=None, wfdb=None, excdb=None): """ Compute the spin matrix elements for excitons. @@ -21,8 +22,8 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut Path to the directory containing the SAVE folder. Default is '.'. bse_dir : str, optional Directory containing the BSE data. Default is 'SAVE'. - iqpt : int, optional - Index of the q-point for which the exciton spin is computed. Default is 1. + iqpt : int or list/array of ints, optional + Index/indices of the q-point(s) for which the exciton spin is computed. Default is 1 (gamma). nstates : int, optional Number of exciton states to consider. If -1, all states are included. Default is -1. sz : array_like, optional @@ -31,12 +32,14 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut 'b','e', 'h'. If 'b' total spin is computed. 'e'/'h' for only electron/hole spin. Default is 'b' (total spin) - if latticedb, wfdb, excitondb are not None, then the provided db is used and is not read again + if latticedb, wfdb, excdb are not None, then the provided db is used and is not read again + excdb is None or single YamboExcitonDB or list of YamboExcitonDB. If this is given, then iqpt is ignored + # # Returns ------- exe_Sz : ndarray - Spin matrix elements for excitons with shape (nstates, nstates). + Spin matrix elements for excitons with shape (niq, nstates, nstates). -------------------------------------------------------------------- Example: @@ -61,54 +64,62 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut # # """ - # Filename for the BSE diagonalization data - filename = 'ndb.BS_diago_Q%d' % (iqpt) + ## Check if it single Q or multiple Q's + if np.isscalar(iqpt): iqpt = [iqpt] + else : iqpt = list(iqpt) # Load the lattice database if not lattice: lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) - # Load the exciton database + ## if not excdb: - excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, + 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_WF=True, neigs=nstates)) + else : + if type(excdb) != list: + excdb = [excdb] # Load the wavefunction database if not wfdb: wfdb = YamboWFDB(path=path, latdb=lattice, - bands_range=[np.min(excdb.table[:, 1]) - 1, - np.max(excdb.table[:, 2])]) - ## sanity check - assert np.min(excdb.table[:, 1]) - 1 == wfdb.min_bnd, \ - "wfdb and exciton db are inconsistant (Bands)" - # - assert np.max(excdb.table[:, 2]) == wfdb.min_bnd + wfdb.nbands, \ - "wfdb and exciton db are inconsistant (Bands)" - + bands_range=[np.min(excdb[0].table[:, 1]) - 1, + np.max(excdb[0].table[:, 2])]) # # Ensure the calculation is valid only for spinor wavefunctions assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" - - # Convert the exciton database to kcv format - excdb.convert_to_kcv() - + # # Compute the spin matrix elements in the BZ elec_sz = wfdb.get_spin_m_e_BZ(s_z=sz) + # + exe_Sz = [] + for ix in range(len(excdb)): + ## sanity check + assert np.min(excdb[ix].table[:, 1]) - 1 == wfdb.min_bnd, \ + "wfdb and exciton db are inconsistant (Bands)" + ## sanity check + assert np.max(excdb[ix].table[:, 2]) == wfdb.min_bnd + wfdb.nbands, \ + "wfdb and exciton db are inconsistant (Bands)" - # Get the exciton q-point in Cartesian coordinates - excQpt = excdb.car_qpoint + # Convert the exciton database to kcv format + excdb[ix].convert_to_kcv() + # Get the exciton q-point in Cartesian coordinates + excQpt = excdb[ix].car_qpoint - # Convert the q-point to crystal coordinates - excQpt = lattice.lat @ excQpt + # Convert the q-point to crystal coordinates + excQpt = lattice.lat @ excQpt - # Compute the exciton spin matrix elements - exe_Sz = exciton_X_matelem(excQpt, np.array([0, 0, 0]), excdb.eigenvectors, - excdb.eigenvectors, elec_sz[None, ...], wfdb.kBZ, - diagonal_only=False,contribution=contribution) + # Compute the exciton spin matrix elements + exe_Sz.append(exciton_X_matelem(excQpt, np.array([0, 0, 0]), excdb[ix].eigenvectors, + excdb[ix].eigenvectors, elec_sz[None, ...], wfdb.kBZ, + diagonal_only=False,contribution=contribution)) # Print a note about the spin matrix #print("Note: This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") - - return exe_Sz + + return np.array(exe_Sz) From b7eb79c5160b45f91eaed03ea82ff825ed86b102 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 6 Feb 2025 23:55:36 +0100 Subject: [PATCH 12/73] dim mismatch in ex spin --- yambopy/bse/exciton_spin.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index a26e7ecf..3c4867ac 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -116,7 +116,7 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut # Compute the exciton spin matrix elements exe_Sz.append(exciton_X_matelem(excQpt, np.array([0, 0, 0]), excdb[ix].eigenvectors, excdb[ix].eigenvectors, elec_sz[None, ...], wfdb.kBZ, - diagonal_only=False,contribution=contribution)) + diagonal_only=False,contribution=contribution)[0]) # Print a note about the spin matrix #print("Note: This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") From 4eeb635bc71f2a08d00b1ef153acffe0777cf82f Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 7 Feb 2025 09:58:37 +0100 Subject: [PATCH 13/73] Improve Exciton spin --- yambopy/bse/exciton_spin.py | 215 ++++++++++++++++++----------- yambopy/dbs/excitondb.py | 14 +- yambopy/dbs/wfdb.py | 35 ++++- yambopy/tools/degeneracy_finder.py | 52 +++++++ 4 files changed, 222 insertions(+), 94 deletions(-) create mode 100644 yambopy/tools/degeneracy_finder.py diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 3c4867ac..25682750 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -4,10 +4,10 @@ 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(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribution='b', - sz=0.5 * np.array([[1, 0], [0, -1]]), - lattice=None, wfdb=None, excdb=None): + +def compute_exciton_spin(lattice, excdb, wfdb, elec_sz, contribution='b',diagonal=False): """ Compute the spin matrix elements for excitons. @@ -18,108 +18,159 @@ def compute_exciton_spin(path='.', bse_dir='SAVE', iqpt=1, nstates=-1, contribut Parameters ---------- - path : str, optional - Path to the directory containing the SAVE folder. Default is '.'. - bse_dir : str, optional - Directory containing the BSE data. Default is 'SAVE'. - iqpt : int or list/array of ints, optional - Index/indices of the q-point(s) for which the exciton spin is computed. Default is 1 (gamma). - nstates : int, optional - Number of exciton states to consider. If -1, all states are included. Default is -1. - sz : array_like, optional - Spin-z operator matrix in the basis of spinor wavefunctions. Default is 0.5 * [[1, 0], [0, -1]]. - contribution : char, optional - 'b','e', 'h'. If 'b' total spin is computed. 'e'/'h' for only electron/hole spin. - Default is 'b' (total spin) - - if latticedb, wfdb, excdb are not None, then the provided db is used and is not read again - excdb is None or single YamboExcitonDB or list of YamboExcitonDB. If this is given, then iqpt is ignored - # - # + 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 (niq, nstates, nstates). - - -------------------------------------------------------------------- - Example: - -------------------------------------------------------------------- - import numpy as np - from yambopy.bse.exciton_spin import compute_exciton_spin - # + Spin matrix elements for excitons with shape (nstates, nstates). + """ # - # Total spin - Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=2) + # Ensure the calculation is valid only for spinor wavefunctions + assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" # - # Only Electron spin - Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=2,contribution='e') + # 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)" # - # Only Hole spin - Sz_exe = compute_exciton_spin(bse_dir='GW_BSE',nstates=2,contribution='h') + 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 # - # The first two states are degenerate so, we diagonalize in 2x2 sub block - w = np.linalg.eigvals(Sz_exe) - print(w) ## spin values of excitons + # Convert the q-point to crystal coordinates + excQpt = lattice.lat @ excQpt # + # Compute the exciton spin matrix elements + exe_Sz = exciton_X_matelem(excQpt, np.array([0, 0, 0]), Akcv, + Akcv, elec_sz[None, ...], wfdb.kBZ, + diagonal_only=diagonal,contribution=contribution) # + return exe_Sz + + + + + +def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, + nstates=-1, contribution='b', degen_tol = 1e-3 + sz=0.5 * np.array([[1, 0], [0, -1]]), + return_dbs_and_spin=True): """ + Compute the spin matrix elements ⟨S'|S_z|S⟩ for excitons. + + This function calculates the spin matrix elements for excitons using + wavefunctions and spin operators. Easy to use interface. Use + compute_exciton_spin() incase you already have those db's + + Parameters + ---------- + path : str, optional + Path to the directory containing the `SAVE` folder. Default is `.`. + bse_dir : str, optional + Directory containing the BSE (Bethe-Salpeter Equation) data. Default is `'SAVE'`. + iqpt : int or list/array of ints, optional + Index or indices of the q-point(s) for which the exciton spin is computed. + Default is `1` (Gamma point). + nstates : int, optional + Number of exciton states to consider. If `-1`, all states are included. Default is `-1`. + contribution : {'b', 'e', 'h'}, optional + Specifies which contribution to compute: + - `'b'`: Total exciton spin (default). + - `'e'`: Electron spin only. + - `'h'`: Hole spin only. + degen_tol : float, optional + Degeneracy tolerance for excitons in eV. Default is `1e-3` eV. + sz : ndarray, optional + Spin-z operator matrix in the basis of spinor wavefunctions. + Default is `0.5 * np.array([[1, 0], [0, -1]])`. + return_dbs_and_spin : bool, optional + return [latticedb, wfdb, excdb (s), elec_spin_matrix] + Default is True + Returns + ------- + exe_Sz : ndarray + Spin matrix elements for excitons with shape `(niq, nstates, nstates)`, + where `niq` is the number of q-points. + if return_dbs_and_spin is True, will also return + [latticedb, wfdb, excdb (s), elec_spin_matrix] + 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 - if not lattice: - lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) - - ## - if not excdb: - 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)) - else : - if type(excdb) != list: - excdb = [excdb] - + 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 - if not wfdb: - wfdb = YamboWFDB(path=path, latdb=lattice, - bands_range=[np.min(excdb[0].table[:, 1]) - 1, - np.max(excdb[0].table[:, 2])]) - # - # Ensure the calculation is valid only for spinor wavefunctions - assert wfdb.nspinor == 2, "Makes sense only for nspinor = 2" + 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 ix in range(len(excdb)): - ## sanity check - assert np.min(excdb[ix].table[:, 1]) - 1 == wfdb.min_bnd, \ - "wfdb and exciton db are inconsistant (Bands)" - ## sanity check - assert np.max(excdb[ix].table[:, 2]) == wfdb.min_bnd + wfdb.nbands, \ - "wfdb and exciton db are inconsistant (Bands)" + for ixdb in excdb: + exe_Sz.append(compute_exciton_spin(lattice, ixdb, + wfdb, elec_sz, + contribution=contribution, + diagonal=False)) + # + exe_Sz = get_spinvals(spin_matrix, eigenvalues, atol=degen_tol) + if return_dbs_and_spin : return np.array(exe_Sz),[lattice, wfdb, excdb, elec_sz] + else : return np.array(exe_Sz) + - # Convert the exciton database to kcv format - excdb[ix].convert_to_kcv() - # Get the exciton q-point in Cartesian coordinates - excQpt = excdb[ix].car_qpoint - # Convert the q-point to crystal coordinates - excQpt = lattice.lat @ excQpt +def get_spinvals(spin_matrix, eigenvalues, atol=1e-3, rtol=1e-3): + degen_idx = find_degeneracy_evs(eigenvalues) + spins = [] + for id in degen_idx: + w = np.linalg.eigvals(spin_matrix[find_degeneracy_evs,:][:,find_degeneracy_evs]) + spins.append(w) + return spins - # Compute the exciton spin matrix elements - exe_Sz.append(exciton_X_matelem(excQpt, np.array([0, 0, 0]), excdb[ix].eigenvectors, - excdb[ix].eigenvectors, elec_sz[None, ...], wfdb.kBZ, - diagonal_only=False,contribution=contribution)[0]) - # Print a note about the spin matrix - #print("Note: This is a spin matrix. Diagonalize the matrix in degenerate subspace to get spin values.") - - return np.array(exe_Sz) diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index f00d1ded..c4bb9224 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -82,9 +82,6 @@ def __init__(self,lattice,Qpt,eigenvalues,l_residual,r_residual,spin_pol='no',ca self.table = table self.eigenvectors = eigenvectors self.spin_pol = spin_pol - self.form = 'bs_tab' ## ('bs_tab or 'kcv') - # if form is bs_tab, them eigenvectors are stored in (nexe,BS_table) - # if form == 'kcv', they are stored as (nexe, nk, nc, nv) @classmethod def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True, neigs=-1): @@ -237,15 +234,17 @@ def write_sorted(self,prefix='yambo'): for i,n in sort_i: f.write("%3d %12.8lf %12.8e\n"%(n+1,eig[n],i)) - def convert_to_kcv(self): + def get_Akcv(self): """ Convert eigenvectors from (neigs,BS_table) -> (neigs,k,c,v) For now, only works for nspin = 1. nspinor = 1/2 also works. """ - if self.form == 'kcv': return - self.form = 'kcv' + assert self.spin_pol == 'no', "Rearrange_Akcv works only for nspin = 1" # + tmp_akcv = getattr(self, 'Akcv', None) + if tmp_akcv is not None: return tmp_akcv + # if self.eigenvectors is None: return None eig_wfcs = self.eigenvectors # @@ -278,7 +277,8 @@ def convert_to_kcv(self): # else : # eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) # - self.eigenvectors = eig_wfcs_returned + self.Akcv = eig_wfcs_returned + return self.Akcv def get_nondegenerate(self,eps=1e-4): """ diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 16741e13..8b61fa28 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -70,7 +70,6 @@ def __init__(self, path=None, save='SAVE', filename='ns.wf', bands_range=[], lat self.path = os.path.join(path, save) self.filename = filename - self.wf_expanded = False ## if true, then wfc's are expanded over full BZ # Read wavefunctions self.read(bands_range=bands_range, latdb=latdb) @@ -445,8 +444,9 @@ def expand_fullBZ(self): ------- None """ - if self.wf_expanded: return - else : self.wf_expanded = True + # + if getattr(self, 'wf_bz', None) is not None: return + # kpt_idx = self.ydb.kpoints_indexes sym_idx = self.ydb.symmetry_indexes nkBZ = len(sym_idx) @@ -487,6 +487,8 @@ def get_BZ_wf(self, ik): Returns: list: Wavefunctions at ik (nspin,nbands,nspinor,ngvec) and G-vectors in crystal coordinates (ngvec,3). """ + if getattr(self, 'wf_bz', None) is None: self.expand_fullBZ() + # return [self.wf_bz[ik][..., :self.ngBZ[ik]], self.g_bz[ik, :self.ngBZ[ik], :]] def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): @@ -512,28 +514,51 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): ## isym >= len(symm_mat)/(1+int(time_rev)) must be timereversal symmetries ## X -> Rx + tau, R matrices are given in symm_mat, tau are frac_vec, time_rev is bool ## (nsym, nk, nspin, Rk_bnd, k_bnd) - if not self.wf_expanded: self.expand_fullBZ() + if getattr(self, 'wf_bz', None) is None: self.expand_fullBZ() + ## + ## Check if already computed for SAVE symetries + dmat_save = getattr(self, 'Dmat', None) + # + # + is_save_symm = False if symm_mat is None or frac_vec is None or time_rev is None: + is_save_symm = True + # + if is_save_symm: + ## if dmats are already computed, return those + if dmat_save is not None : return dmat_save + # symm_mat = self.ydb.sym_car frac_vec = np.zeros((len(symm_mat),3),dtype=symm_mat.dtype) time_rev = int(np.rint(self.ydb.time_rev)) - + # ktree = build_ktree(self.kBZ) Dmat = [] nsym = len(symm_mat) + # assert nsym == len(frac_vec), "The number for frac translation must be same as Rotation matrices" for ik in tqdm(range(self.nkBZ), desc="Dmat"): + # wfc_k, gvec_k = self.get_BZ_wf(ik) kvec = self.get_BZ_kpt(ik) + # for isym in range(nsym): trev = (isym >= nsym/(1+int(time_rev))) + # ## Compute U(R)\psi_k + # Rk, wfc_Rk, gvec_Rk = self.apply_symm(kvec, wfc_k, gvec_k, trev, symm_mat[isym], frac_vec[isym]) idx = find_kpt(ktree, Rk) + # ## get Rk wfc stored + # w_rk, g_rk = self.get_BZ_wf(idx) Dmat.append(wfc_inner_product(self.get_BZ_kpt(idx),w_rk, g_rk, Rk, wfc_Rk, gvec_Rk)) + # Dmat = np.array(Dmat).reshape(self.nkBZ, nsym, self.nspin, self.nbands, self.nbands).transpose(1,0,2,3,4) + # + if is_save_symm: self.Dmat = Dmat + # return Dmat def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gtree=-1): diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py new file mode 100644 index 00000000..dc32ce55 --- /dev/null +++ b/yambopy/tools/degeneracy_finder.py @@ -0,0 +1,52 @@ +import numpy as np + +def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): + """ + Identify sets of degenerate eigenstates based on eigenvalues. + + Parameters + ---------- + eigenvalues : array-like + A sorted list or array of eigenvalues in increasing order. + atol : float, optional + Absolute tolerance for degeneracy comparison. Default is 1e-3. + rtol : float, optional + Relative tolerance for degeneracy comparison. Default is 1e-3. + + Returns + ------- + list of lists + A list where each sublist contains the indices of degenerate eigenstates. + + Raises + ------ + ValueError + If `eigenvalues` is empty or not a valid array. + """ + # Input validation + eigenvalues = np.asarray(eigenvalues) + if eigenvalues.size == 0: + raise ValueError("Input eigenvalues must not be empty.") + if atol < 0 or rtol < 0: + raise ValueError("Tolerances `atol` and `rtol` must be non-negative.") + + # Compute differences between consecutive eigenvalues + diffs = np.diff(eigenvalues) + tolerance = atol + rtol * np.abs(eigenvalues[:-1]) + + # Identify where the differences exceed the tolerance + split_indices = np.where(diffs > tolerance)[0] + + # Group indices of degenerate states + degen_sets = np.split(np.arange(len(eigenvalues)), split_indices + 1) + + # Convert numpy arrays to lists for consistency + degen_sets = [group for group in degen_sets] + + return degen_sets + + +if __name__ == '__main__': + eigenvalues = [1.0, 1.001, 1.002, 2.0, 2.001, 3.0] + degen_sets = find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3) + print(degen_sets) From 8df3f534b6ec74d999738f3be36c4186e7f84669 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 7 Feb 2025 10:08:06 +0100 Subject: [PATCH 14/73] Fix error --- yambopy/bse/exciton_spin.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 25682750..54480eb4 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -72,7 +72,7 @@ def compute_exciton_spin(lattice, excdb, wfdb, elec_sz, contribution='b',diagona def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, - nstates=-1, contribution='b', degen_tol = 1e-3 + nstates=-1, contribution='b', degen_tol = 1e-3, sz=0.5 * np.array([[1, 0], [0, -1]]), return_dbs_and_spin=True): """ From f5d6623bad0c88d4aa181d99d5dc9ef36d3e77a3 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 7 Feb 2025 10:27:28 +0100 Subject: [PATCH 15/73] Fix errors in exciton spin --- yambopy/bse/exciton_spin.py | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 54480eb4..a5428556 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -65,7 +65,7 @@ def compute_exciton_spin(lattice, excdb, wfdb, elec_sz, contribution='b',diagona Akcv, elec_sz[None, ...], wfdb.kBZ, diagonal_only=diagonal,contribution=contribution) # - return exe_Sz + return exe_Sz[0] @@ -153,14 +153,18 @@ def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, # exe_Sz = [] for ixdb in excdb: - exe_Sz.append(compute_exciton_spin(lattice, ixdb, + smat = compute_exciton_spin(lattice, ixdb, wfdb, elec_sz, contribution=contribution, - diagonal=False)) + 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 = get_spinvals(spin_matrix, eigenvalues, atol=degen_tol) - if return_dbs_and_spin : return np.array(exe_Sz),[lattice, wfdb, excdb, elec_sz] - else : return np.array(exe_Sz) + exe_Sz = np.array(exe_Sz) + if return_dbs_and_spin : return exe_Sz,[lattice, wfdb, excdb, elec_sz] + else : return exe_Sz @@ -168,7 +172,7 @@ def get_spinvals(spin_matrix, eigenvalues, atol=1e-3, rtol=1e-3): degen_idx = find_degeneracy_evs(eigenvalues) spins = [] for id in degen_idx: - w = np.linalg.eigvals(spin_matrix[find_degeneracy_evs,:][:,find_degeneracy_evs]) + w = np.linalg.eigvals(spin_matrix[id,:][:,id]) spins.append(w) return spins From bf6501de6344365f554838f83908ff346fe7abd8 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 7 Feb 2025 16:51:18 +0100 Subject: [PATCH 16/73] Add some functions for lazy loading in letz --- yambopy/letzelphc_interface/lelph2y.py | 15 ++-- yambopy/letzelphc_interface/lelphcdb.py | 93 ++++++++++++++++++++++++- 2 files changed, 102 insertions(+), 6 deletions(-) diff --git a/yambopy/letzelphc_interface/lelph2y.py b/yambopy/letzelphc_interface/lelph2y.py index 5657231e..37c7595c 100644 --- a/yambopy/letzelphc_interface/lelph2y.py +++ b/yambopy/letzelphc_interface/lelph2y.py @@ -50,9 +50,14 @@ def __init__(self,OBJ,code,SAVE_path,OUT_path=None): self.get_yambo_header_variables(SAVE_path) # Get el-ph data from external code - match code: - case 'lelphc': self.get_elph_variables_LELPHC(OBJ) - case _: raise NotImplementedError("Code %s not found or implemented"%code) + # + # NM : using match will enforce python 3.10. some HPC's still use <= 3.8 + # Commenting it out until it gets old enough. + # match code: + # case 'lelphc': self.get_elph_variables_LELPHC(OBJ) + # case _: raise NotImplementedError("Code %s not found or implemented"%code) + if code.strip() == 'lelphc': self.get_elph_variables_LELPHC(OBJ) + else: raise NotImplementedError("Code %s not found or implemented"%code) if not os.path.isdir(OUT_path): os.mkdir(OUT_path) @@ -137,7 +142,7 @@ def write_header(self,OUT_path): dbs.createDimension('D_%010d'%1,1) dbs.createDimension('D_%010d'%2,2) dbs.createDimension('D_%010d'%4,4) - for value in [self.natoms,self.nkpoints_ibz,len_pars,self.nqpoints_bz]: + for value in [self.natoms,self.nkpoints_ibz,len_pars,self.nqpoints_bz, self.nkpoints_bz]: if value not in [1,2,3,4]: try: dbs.createDimension('D_%010d'%value,value) except RuntimeError: pass # This is when one of the dimensions already exists @@ -154,7 +159,7 @@ def write_header(self,OUT_path): PARS = dbs.createVariable('PARS',netcdftype('f',p), ('D_%.10d'%len_pars)) dbs.createVariable('MAX_PH_FREQ',netcdftype('f',p), ('D_%.10d'%1)) dbs.createVariable('PH_Q',netcdftype('f',p), ('D_%.10d'%3, 'D_%.10d'%self.nqpoints_bz)) - dbs.createVariable('PH_K',netcdftype('f',p), ('D_%.10d'%3, 'D_%.10d'%self.nqpoints_bz)) + dbs.createVariable('PH_K',netcdftype('f',p), ('D_%.10d'%3, 'D_%.10d'%self.nkpoints_bz)) # Fill variables dbs['HEAD_VERSION'][:] = self.head_version diff --git a/yambopy/letzelphc_interface/lelphcdb.py b/yambopy/letzelphc_interface/lelphcdb.py index e25c1c7d..f75d81bf 100644 --- a/yambopy/letzelphc_interface/lelphcdb.py +++ b/yambopy/letzelphc_interface/lelphcdb.py @@ -2,6 +2,7 @@ from netCDF4 import Dataset from yambopy.tools.string import marquee from yambopy.units import ha2ev +from yambopy.kpoints import build_ktree, find_kpt class LetzElphElectronPhononDB(): """ @@ -39,6 +40,7 @@ def __init__(self,filename,read_all=True,div_by_energies=True): try: database = Dataset(filename) except: raise FileNotFoundError("error opening %s in LetzElphElectronPhononDB"%filename) + self.filename = filename # Read DB dimensions self.nb1 = database.dimensions['initial_band'].size self.nb2 = database.dimensions['final_band_PH_abs'].size @@ -49,10 +51,23 @@ def __init__(self,filename,read_all=True,div_by_energies=True): self.ns = database.dimensions['nspin'].size self.nsym = database.dimensions['nsym_ph'].size + conv = database['convention'][...].data + convlist = [iconv.decode('utf-8') for iconv in conv] + conv = '' + for iconv in convlist: conv = conv + iconv + stard_conv = False + if conv.strip() == 'standard': + print("Convention used in Letzelphc : k -> k+q (standard)") + else: + print("Convention used in Letzelphc : k-q -> k (yambo)") + self.convention = conv.strip() + # # Read DB self.kpoints = database.variables['kpoints'][:] self.qpoints = database.variables['qpoints'][:] self.bands = database.variables['bands'][:] + self.ktree = build_ktree(self.kpoints) + self.qtree = build_ktree(self.qpoints) self.ph_energies = database.variables['FREQ'][:]*(ha2ev/2.) # Energy units are in Rydberg self.check_energies() @@ -128,6 +143,82 @@ def scale_g(self,dvscf): g[iq,:,inu,:,:,:] = dvscf[iq,:,inu,:,:,:]/np.sqrt(2.*ph_E) return dvscf + def read_read_iq(self,iq, bands_range=[], database=None, convention='yambo'): + """ + Reads the electron-phonon matrix elements and phonon eigenvectors for a specific q-point index. + + If the data is already loaded in memory, it returns the corresponding array slice. Otherwise, + it reads from the database without storing the data in memory. + + This function reads data for a single q-point instead of the entire dataset, which is useful + for handling large databases efficiently. + + Parameters + ---------- + iq : int + Index of the q-point. + bands_range : list, optional + Specifies the range of bands to read. The start index follows Python indexing (starting from 0), + and the end index is excluded. If not provided, it defaults to the minimum and maximum bands available. + database : Dataset, optional + If provided, the function will use this open dataset instead of opening the file again. + convention : str, optional + Defines the convention used for electron-phonon matrix elements. + - 'yambo': Outputs \. + - Any other value: Outputs \. + + Returns + ------- + tuple + A tuple containing: + - ph_eigenvectors : ndarray + The phonon eigenvectors. + - ph_elph_me : ndarray + The electron-phonon matrix elements with the specified convention. + """ + # + if len(bands_range) == 0: + bands_range = [min(self.bands)-1,max(self.bands)] + min_bnd = min(bands_range) + max_bnd = max(bands_range) + nbnds = max_bnd - min_bnd + assert (min_bnd >= min(self.bands)-1) + assert (max_bnd <= max(self.bands)) + start_bnd_idx = 1+min_bnd - min(self.bands) + end_bnd = start_bnd_idx + nbnds + + # self.ph_eigenvectors , self.gkkp + if hasattr(self, 'ph_eigenvectors'): + ph_eigs = self.ph_eigenvectors[iq] + eph_mat = self.gkkp[iq, :, :, :, start_bnd_idx:end_bnd, start_bnd_idx:end_bnd ] + else : + ## else we load from the file + close_file = False + if not database : + close_file = True + database = Dataset(filename,'r') + eph_mat = database['elph_mat'][iq, :, :, :, start_bnd_idx:end_bnd, start_bnd_idx:end_bnd, :].data + ph_eigs = database['POLARIZATION_VECTORS'][iq,...].data + eph_mat = eph_mat[...,0] + 1j*eph_mat[...,1] + ph_eigs = ph_eigs[...,0] + 1j*ph_eigs[...,1] + if close_file :database.close() + return [ph_eigs, self.change_convention(self.qpoints[iq],eph_mat, convention)] + + def change_convention(self, qpt, elph_iq, convention='yambo'): + """ + ## (nk,....) + qpt in crystal coordinates + if convention == 'yambo', it will give '' + else + # This returns a view and not a copy + """ + if convention.strip() != 'yambo': convention = 'standard' + if self.convention == convention: return elph_iq + if convention == 'standard': factor = 1.0 + else factor = -1.0 + idx_q = find_kindx(self.ktree, factor*qpt[None, :] + self.kpoints) + return elph_iq[idx_q, ...] + def __str__(self,verbose=False): lines = []; app = lines.append @@ -155,4 +246,4 @@ def __str__(self,verbose=False): app('Eigenvectors and matrix elements not loaded') app('-----------------------------------') - return "\n".join(lines) \ No newline at end of file + return "\n".join(lines) From 9db1a1108d9322c65068e1f1a62a0437e557f743 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 7 Feb 2025 17:04:31 +0100 Subject: [PATCH 17/73] Lazy loading for letz --- yambopy/letzelphc_interface/lelphcdb.py | 30 +++++++++++++++++-------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/yambopy/letzelphc_interface/lelphcdb.py b/yambopy/letzelphc_interface/lelphcdb.py index f75d81bf..53b1a7f0 100644 --- a/yambopy/letzelphc_interface/lelphcdb.py +++ b/yambopy/letzelphc_interface/lelphcdb.py @@ -143,7 +143,7 @@ def scale_g(self,dvscf): g[iq,:,inu,:,:,:] = dvscf[iq,:,inu,:,:,:]/np.sqrt(2.*ph_E) return dvscf - def read_read_iq(self,iq, bands_range=[], database=None, convention='yambo'): + def read_iq(self,iq, bands_range=[], database=None, convention='yambo'): """ Reads the electron-phonon matrix elements and phonon eigenvectors for a specific q-point index. @@ -196,7 +196,7 @@ def read_read_iq(self,iq, bands_range=[], database=None, convention='yambo'): close_file = False if not database : close_file = True - database = Dataset(filename,'r') + database = Dataset(self.filename,'r') eph_mat = database['elph_mat'][iq, :, :, :, start_bnd_idx:end_bnd, start_bnd_idx:end_bnd, :].data ph_eigs = database['POLARIZATION_VECTORS'][iq,...].data eph_mat = eph_mat[...,0] + 1j*eph_mat[...,1] @@ -206,17 +206,29 @@ def read_read_iq(self,iq, bands_range=[], database=None, convention='yambo'): def change_convention(self, qpt, elph_iq, convention='yambo'): """ - ## (nk,....) - qpt in crystal coordinates - if convention == 'yambo', it will give '' - else - # This returns a view and not a copy + Adjusts the convention of the electron-phonon matrix elements. + + Parameters + ---------- + qpt : ndarray + The q-point in crystal coordinates. + elph_iq : ndarray + The electron-phonon matrix elements. + convention : str, optional + Defines the output format: + - 'yambo': Outputs \. + - Any other value: Outputs \. + + Returns + ------- + ndarray + The electron-phonon matrix elements in the desired convention. The returned array is a view, not a copy. """ if convention.strip() != 'yambo': convention = 'standard' if self.convention == convention: return elph_iq if convention == 'standard': factor = 1.0 - else factor = -1.0 - idx_q = find_kindx(self.ktree, factor*qpt[None, :] + self.kpoints) + else :factor = -1.0 + idx_q = find_kpt(self.ktree, factor*qpt[None, :] + self.kpoints) return elph_iq[idx_q, ...] def __str__(self,verbose=False): From 2dfc557c2342a15a5d740d92889ffec9359d6a1e Mon Sep 17 00:00:00 2001 From: Murali Date: Sat, 8 Feb 2025 10:07:07 +0100 Subject: [PATCH 18/73] Dmat can be computed without wf expansion --- yambopy/dbs/wfdb.py | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 8b61fa28..225d98c0 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -514,7 +514,8 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): ## isym >= len(symm_mat)/(1+int(time_rev)) must be timereversal symmetries ## X -> Rx + tau, R matrices are given in symm_mat, tau are frac_vec, time_rev is bool ## (nsym, nk, nspin, Rk_bnd, k_bnd) - if getattr(self, 'wf_bz', None) is None: self.expand_fullBZ() + expand_wf_present = True + if getattr(self, 'wf_bz', None) is None: expand_wf_present = False ## ## Check if already computed for SAVE symetries dmat_save = getattr(self, 'Dmat', None) @@ -535,13 +536,21 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): ktree = build_ktree(self.kBZ) Dmat = [] nsym = len(symm_mat) + kpt_idx = self.ydb.kpoints_indexes + sym_idx = self.ydb.symmetry_indexes # assert nsym == len(frac_vec), "The number for frac translation must be same as Rotation matrices" for ik in tqdm(range(self.nkBZ), desc="Dmat"): - # - wfc_k, gvec_k = self.get_BZ_wf(ik) - kvec = self.get_BZ_kpt(ik) - # + # IN case the wfc are already expanded, load them + if expand_wf_present: + wfc_k, gvec_k = self.get_BZ_wf(ik) + kvec = self.get_BZ_kpt(ik) + else: + ## else rotate them + ixk = kpt_idx[ik] + isk = sym_idx[ik] + kvec, wfc_k, gvec_k = self.rotate_wfc(ixk, isk) + for isym in range(nsym): trev = (isym >= nsym/(1+int(time_rev))) # @@ -550,10 +559,16 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): Rk, wfc_Rk, gvec_Rk = self.apply_symm(kvec, wfc_k, gvec_k, trev, symm_mat[isym], frac_vec[isym]) idx = find_kpt(ktree, Rk) # - ## get Rk wfc stored - # - w_rk, g_rk = self.get_BZ_wf(idx) - Dmat.append(wfc_inner_product(self.get_BZ_kpt(idx),w_rk, g_rk, Rk, wfc_Rk, gvec_Rk)) + ## get Rk wfc + # in case stored + if expand_wf_present: + w_rk, g_rk = self.get_BZ_wf(idx) + k_rk = self.get_BZ_kpt(idx) + else : + iktmp = kpt_idx[idx] + istmp = sym_idx[idx] + k_rk, w_rk, g_rk = self.rotate_wfc(iktmp, istmp) + Dmat.append(wfc_inner_product(k_rk, w_rk, g_rk, Rk, wfc_Rk, gvec_Rk)) # Dmat = np.array(Dmat).reshape(self.nkBZ, nsym, self.nspin, self.nbands, self.nbands).transpose(1,0,2,3,4) # From 734ba69b2544174af3de6c249de88beb3b1db76d Mon Sep 17 00:00:00 2001 From: Murali Date: Sat, 8 Feb 2025 10:28:45 +0100 Subject: [PATCH 19/73] Fix messed up Dmat var --- yambopy/dbs/wfdb.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 225d98c0..2a7bdfad 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -518,7 +518,7 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): if getattr(self, 'wf_bz', None) is None: expand_wf_present = False ## ## Check if already computed for SAVE symetries - dmat_save = getattr(self, 'Dmat', None) + dmat_save = getattr(self, 'save_Dmat', None) # # is_save_symm = False @@ -572,7 +572,7 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): # Dmat = np.array(Dmat).reshape(self.nkBZ, nsym, self.nspin, self.nbands, self.nbands).transpose(1,0,2,3,4) # - if is_save_symm: self.Dmat = Dmat + if is_save_symm: self.save_Dmat = Dmat # return Dmat From 78cfe97cad2e5d7226a9fa8e8ae94cc616c516ec Mon Sep 17 00:00:00 2001 From: Murali Date: Sat, 8 Feb 2025 10:36:33 +0100 Subject: [PATCH 20/73] Fix error when merging bux fixes branch --- yambopy/letzelphc_interface/lelph2y.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/letzelphc_interface/lelph2y.py b/yambopy/letzelphc_interface/lelph2y.py index 42a20d10..126b6ee5 100644 --- a/yambopy/letzelphc_interface/lelph2y.py +++ b/yambopy/letzelphc_interface/lelph2y.py @@ -143,7 +143,7 @@ def write_header(self,OUT_path): dbs.createDimension('D_%010d'%2,2) dbs.createDimension('D_%010d'%4,4) # - for value in [self.natoms,self.nkpoints_ibz,len_pars,self.nqpoints_bz,self.nk_points_bz]: + for value in [self.natoms,self.nkpoints_ibz,len_pars,self.nqpoints_bz,self.nkpoints_bz]: if value not in [1,2,3,4]: try: dbs.createDimension('D_%010d'%value,value) except RuntimeError: pass # This is when one of the dimensions already exists From 06ed17791b958b989d1c5213db11b6cd43fe0ce9 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 9 Feb 2025 21:04:40 +0100 Subject: [PATCH 21/73] reduce python dependency version --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 732a9f93..0d52f76e 100644 --- a/pyproject.toml +++ b/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+)", From 114da3d9564299e276659ead29d388e246e980d0 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 9 Feb 2025 22:37:28 +0100 Subject: [PATCH 22/73] Add comments --- yambopy/bse/exciton_matrix_elements.py | 3 ++- yambopy/letzelphc_interface/lelphcdb.py | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index fb0be3b8..6364282d 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -28,7 +28,8 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di Ak : array_like Wavefunction coefficients for k (ket wfc) with shape (n_exe_states, nk, nc, nv). Omn : array_like - Matrix elements of the operator O in the basis of electronic states with shape (nlambda, nk, nm, nn). + Matrix elements of the operator O in the basis of electronic states with shape (nlambda, nk, m_bnd, n_bnd). + ie Omn = < k+q, m | O(q) | n, k>, where m_bnd and n_bnd are initial and final bands. kpts : array_like K-points used to construct the BSE with shape (nk, 3). contribution : str, optional diff --git a/yambopy/letzelphc_interface/lelphcdb.py b/yambopy/letzelphc_interface/lelphcdb.py index 53b1a7f0..08f642b8 100644 --- a/yambopy/letzelphc_interface/lelphcdb.py +++ b/yambopy/letzelphc_interface/lelphcdb.py @@ -198,6 +198,7 @@ def read_iq(self,iq, bands_range=[], database=None, convention='yambo'): close_file = True database = Dataset(self.filename,'r') eph_mat = database['elph_mat'][iq, :, :, :, start_bnd_idx:end_bnd, start_bnd_idx:end_bnd, :].data + # ( nk, nm, nspin, initial bnd, final bnd) ph_eigs = database['POLARIZATION_VECTORS'][iq,...].data eph_mat = eph_mat[...,0] + 1j*eph_mat[...,1] ph_eigs = ph_eigs[...,0] + 1j*ph_eigs[...,1] From 39b3479e736ee981bc58cd3c5f63a27915f9d112 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 10 Feb 2025 10:06:34 +0100 Subject: [PATCH 23/73] Add more comments --- yambopy/dbs/wfdb.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 245dab81..3364a2e8 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -476,6 +476,8 @@ def expand_fullBZ(self): self.wf_bz = np.zeros((nkBZ, self.nspin, self.nbands, self.nspinor, self.ng),dtype=self.wf.dtype) self.g_bz = 2147483646 + np.zeros((nkBZ,self.ng,3),dtype=int) self.kBZ = np.zeros((nkBZ,3)) + # NM : The reason we want to replace the existing kBZ variable is to make sure we have + # correct rotated kpoint (here they should not differ by a G vector !) self.ngBZ = np.zeros(nkBZ,dtype=int) for i in tqdm(range(nkBZ), desc="Expanding Wavefunctions full BZ"): ik = kpt_idx[i] @@ -487,6 +489,7 @@ def expand_fullBZ(self): self.kBZ[i] = kbz self.wf_bz[i][...,:ng_t] = w_t self.g_bz[i][:ng_t,:] = g_t + self.kBZ_wf = self.kBZ ## Saving for sanity purposes. because these kvecs should not differ by G def get_BZ_kpt(self, ik): """ @@ -498,7 +501,8 @@ def get_BZ_kpt(self, ik): Returns: numpy.ndarray: K-point in crystal coordinates. """ - return self.kBZ[ik] + kpts_tmp = getattr(self, 'kBZ_wf', self.kBZ) + return kpts_tmp[ik] def get_BZ_wf(self, ik): """ @@ -564,6 +568,8 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): assert nsym == len(frac_vec), "The number for frac translation must be same as Rotation matrices" for ik in tqdm(range(self.nkBZ), desc="Dmat"): # IN case the wfc are already expanded, load them + # NM : Be very carefull when dealing with kvectors here, use the correct kvec + # for the wfc (they should not differ by G vector !) if expand_wf_present: wfc_k, gvec_k = self.get_BZ_wf(ik) kvec = self.get_BZ_kpt(ik) From 0a5e6538ce259f2329119688c2e48224e74e965f Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 10 Feb 2025 10:08:08 +0100 Subject: [PATCH 24/73] Add more comments --- yambopy/dbs/wfdb.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 3364a2e8..793c7d8f 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -489,7 +489,6 @@ def expand_fullBZ(self): self.kBZ[i] = kbz self.wf_bz[i][...,:ng_t] = w_t self.g_bz[i][:ng_t,:] = g_t - self.kBZ_wf = self.kBZ ## Saving for sanity purposes. because these kvecs should not differ by G def get_BZ_kpt(self, ik): """ @@ -501,8 +500,7 @@ def get_BZ_kpt(self, ik): Returns: numpy.ndarray: K-point in crystal coordinates. """ - kpts_tmp = getattr(self, 'kBZ_wf', self.kBZ) - return kpts_tmp[ik] + return self.kBZ[ik] def get_BZ_wf(self, ik): """ @@ -568,8 +566,6 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): assert nsym == len(frac_vec), "The number for frac translation must be same as Rotation matrices" for ik in tqdm(range(self.nkBZ), desc="Dmat"): # IN case the wfc are already expanded, load them - # NM : Be very carefull when dealing with kvectors here, use the correct kvec - # for the wfc (they should not differ by G vector !) if expand_wf_present: wfc_k, gvec_k = self.get_BZ_wf(ik) kvec = self.get_BZ_kpt(ik) From 6653649ac575420c5d82d09590e107f654794806 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 10 Feb 2025 10:51:43 +0100 Subject: [PATCH 25/73] fix error in comment --- yambopy/bse/exciton_matrix_elements.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index 6364282d..c334ca65 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -29,7 +29,7 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di Wavefunction coefficients for k (ket wfc) with shape (n_exe_states, nk, nc, nv). Omn : array_like Matrix elements of the operator O in the basis of electronic states with shape (nlambda, nk, m_bnd, n_bnd). - ie Omn = < k+q, m | O(q) | n, k>, where m_bnd and n_bnd are initial and final bands. + ie Omn = < k+q, m | O(q) | n, k>, where m_bnd and n_bnd are final and initial bands, respectively. kpts : array_like K-points used to construct the BSE with shape (nk, 3). contribution : str, optional From 4658c4beb205bb0730c91ab6ef901c22a339aeae Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 11 Feb 2025 18:18:40 +0100 Subject: [PATCH 26/73] Improve degeneracy finder --- yambopy/tools/degeneracy_finder.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py index dc32ce55..71b45ded 100644 --- a/yambopy/tools/degeneracy_finder.py +++ b/yambopy/tools/degeneracy_finder.py @@ -7,7 +7,7 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): Parameters ---------- eigenvalues : array-like - A sorted list or array of eigenvalues in increasing order. + A list or array of eigenvalues atol : float, optional Absolute tolerance for degeneracy comparison. Default is 1e-3. rtol : float, optional @@ -25,6 +25,8 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): """ # Input validation eigenvalues = np.asarray(eigenvalues) + idx_sorted = np.argsort(eigenvalues) + eigenvalues = eigenvalues[idx_sorted] if eigenvalues.size == 0: raise ValueError("Input eigenvalues must not be empty.") if atol < 0 or rtol < 0: @@ -41,12 +43,12 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): degen_sets = np.split(np.arange(len(eigenvalues)), split_indices + 1) # Convert numpy arrays to lists for consistency - degen_sets = [group for group in degen_sets] + degen_sets = [idx_sorted[group] for group in degen_sets] return degen_sets if __name__ == '__main__': - eigenvalues = [1.0, 1.001, 1.002, 2.0, 2.001, 3.0] + eigenvalues = [3.0001, 2.99976, 1.99999, 1.0, 1.001, 1.002, 2.0, 2.001, 3.0] degen_sets = find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3) print(degen_sets) From 6e07308aba883e5cd282b87bf33c7cceb682c144 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 11 Feb 2025 19:37:16 +0100 Subject: [PATCH 27/73] small change in degeneracy finder --- yambopy/tools/degeneracy_finder.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py index 71b45ded..7948401d 100644 --- a/yambopy/tools/degeneracy_finder.py +++ b/yambopy/tools/degeneracy_finder.py @@ -23,15 +23,16 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): ValueError If `eigenvalues` is empty or not a valid array. """ - # Input validation eigenvalues = np.asarray(eigenvalues) - idx_sorted = np.argsort(eigenvalues) - eigenvalues = eigenvalues[idx_sorted] + # if eigenvalues.size == 0: raise ValueError("Input eigenvalues must not be empty.") if atol < 0 or rtol < 0: raise ValueError("Tolerances `atol` and `rtol` must be non-negative.") - + # + idx_sorted = np.argsort(eigenvalues) + eigenvalues = eigenvalues[idx_sorted] + # # Compute differences between consecutive eigenvalues diffs = np.diff(eigenvalues) tolerance = atol + rtol * np.abs(eigenvalues[:-1]) From 3c3c8157f46c325bc7cdfacdef2db3cee3775fd0 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 16 Feb 2025 23:52:53 +0100 Subject: [PATCH 28/73] Avoid some suitations in degeneracy finder. --- yambopy/tools/degeneracy_finder.py | 50 +++++++++++++++++++++++------- 1 file changed, 38 insertions(+), 12 deletions(-) diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py index 7948401d..8b84dcc9 100644 --- a/yambopy/tools/degeneracy_finder.py +++ b/yambopy/tools/degeneracy_finder.py @@ -1,4 +1,4 @@ -import numpy as np +import numpy as np def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): """ @@ -24,29 +24,55 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): If `eigenvalues` is empty or not a valid array. """ eigenvalues = np.asarray(eigenvalues) - # + if eigenvalues.size == 0: raise ValueError("Input eigenvalues must not be empty.") if atol < 0 or rtol < 0: raise ValueError("Tolerances `atol` and `rtol` must be non-negative.") - # + + # Sort eigenvalues and get sorted indices idx_sorted = np.argsort(eigenvalues) eigenvalues = eigenvalues[idx_sorted] - # + # Compute differences between consecutive eigenvalues diffs = np.diff(eigenvalues) tolerance = atol + rtol * np.abs(eigenvalues[:-1]) - + # Identify where the differences exceed the tolerance split_indices = np.where(diffs > tolerance)[0] - + # Group indices of degenerate states - degen_sets = np.split(np.arange(len(eigenvalues)), split_indices + 1) - - # Convert numpy arrays to lists for consistency - degen_sets = [idx_sorted[group] for group in degen_sets] - - return degen_sets + degen_sets = np.split(idx_sorted, split_indices + 1) + + # Further split groups based on the mean of the group + # NM : This is to ensure that the list is a Arithmetic progression + # with d < tol + final_degen_sets = [] + for group in degen_sets: + if len(group) == 0: + continue + group_eigenvalues = eigenvalues[group] + current_group = [group[0]] + current_mean = group_eigenvalues[0] + + for i in range(1, len(group)): + diff = np.abs(group_eigenvalues[i] - current_mean) + tolerance = atol + rtol * np.abs(current_mean) + + if diff <= tolerance: + current_group.append(group[i]) + # Update the mean of the current group incrementally + current_mean = (current_mean * (len(current_group) - 1) + group_eigenvalues[i]) / len(current_group) + else: + final_degen_sets.append(current_group) + current_group = [group[i]] + current_mean = group_eigenvalues[i] + + # Append the last group + if current_group: + final_degen_sets.append(current_group) + + return final_degen_sets if __name__ == '__main__': From 7533570fc2c80c29eb3c7f3e3d80c84f0ff13455 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 17 Feb 2025 00:05:43 +0100 Subject: [PATCH 29/73] fix bug in denegeray finder --- yambopy/tools/degeneracy_finder.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py index 8b84dcc9..f2b42f07 100644 --- a/yambopy/tools/degeneracy_finder.py +++ b/yambopy/tools/degeneracy_finder.py @@ -32,23 +32,24 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): # Sort eigenvalues and get sorted indices idx_sorted = np.argsort(eigenvalues) - eigenvalues = eigenvalues[idx_sorted] + eigenvalues_sorted = eigenvalues[idx_sorted] # Compute differences between consecutive eigenvalues - diffs = np.diff(eigenvalues) - tolerance = atol + rtol * np.abs(eigenvalues[:-1]) + diffs = np.diff(eigenvalues_sorted) + tolerance = atol + rtol * np.abs(eigenvalues_sorted[:-1]) # Identify where the differences exceed the tolerance split_indices = np.where(diffs > tolerance)[0] - # Group indices of degenerate states - degen_sets = np.split(idx_sorted, split_indices + 1) + # Group indices of degenerate states (in sorted order) + degen_sets_sorted = np.split(idx_sorted, split_indices + 1) # Further split groups based on the mean of the group - # NM : This is to ensure that the list is a Arithmetic progression - # with d < tol + # NM : This is to ensure that in case the sequecence if + # Arthematic progession with d < tol, we make sure that the + # values are within its mean final_degen_sets = [] - for group in degen_sets: + for group in degen_sets_sorted: if len(group) == 0: continue group_eigenvalues = eigenvalues[group] @@ -64,13 +65,13 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): # Update the mean of the current group incrementally current_mean = (current_mean * (len(current_group) - 1) + group_eigenvalues[i]) / len(current_group) else: - final_degen_sets.append(current_group) + final_degen_sets.append(np.array(current_group)) current_group = [group[i]] current_mean = group_eigenvalues[i] # Append the last group if current_group: - final_degen_sets.append(current_group) + final_degen_sets.append(np.array(current_group)) return final_degen_sets From 691c5af60353542fbb5d90ea8be95cb433e8850b Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 25 Feb 2025 13:08:35 +0100 Subject: [PATCH 30/73] fix tiny bug in wfdb --- yambopy/dbs/wfdb.py | 1 + 1 file changed, 1 insertion(+) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 793c7d8f..239debdc 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -412,6 +412,7 @@ def apply_symm(self, kvec, wfc_k, gvecs_k, time_rev, sym_mat, frac_vec=np.array( su_mat = su2_mat(sym_mat, time_rev).astype(wfc_k.dtype) #'ij,sbjg->sbig', su_mat, wfc_k wfc_rot = su_mat[None,None,:,:]@wfc_k + else : wfc_rot = wfc_k.copy() # Add phase due to fractional translation Rkvec = sym_red.T @ kvec From a01600108a62441b4f4c942d9ec5938425d121f9 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 25 Mar 2025 17:17:09 +0100 Subject: [PATCH 31/73] nonTDA check symm --- yambopy/bse/exciton_matrix_elements.py | 1 + yambopy/bse/rotate_excitonwf.py | 18 +++++++++++++++--- yambopy/dbs/excitondb.py | 15 ++++++++------- 3 files changed, 24 insertions(+), 10 deletions(-) diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index c334ca65..38cb7ffd 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -51,6 +51,7 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di # Number of arbitrary parameters (lambda) in the Omn matrix nlambda = Omn.shape[0] # + assert len(Akq.shape) == 4, "Works only with TDA." # Shape of the wavefunction coefficients n_exe_states, nk, nc, nv = Akq.shape # diff --git a/yambopy/bse/rotate_excitonwf.py b/yambopy/bse/rotate_excitonwf.py index 0e9ca918..72a2ad8e 100644 --- a/yambopy/bse/rotate_excitonwf.py +++ b/yambopy/bse/rotate_excitonwf.py @@ -14,6 +14,7 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non ---------- Ak : array_like Exciton wavefunction coefficients with shape (n_exe_states, nk, nc, nv). + In case of non-TDA, it is (n_exe_states, 2, nk, nc, nv) symm_mat_red : array_like Symmetry matrix in reduced coordinates with shape (3, 3). kpoints : array_like @@ -34,8 +35,11 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non """ # Initialize the rotated Ak array rot_Ak = np.zeros(Ak.shape, dtype=Ak.dtype) - nk, nc, nv = Ak.shape[1:] + # Check TDA + tda = True + if len(Ak.shape) != 4 : tda = False + nk, nc, nv = Ak.shape[-3:] # Build a k-point tree if not provided if not ktree: ktree = build_ktree(kpoints) @@ -56,7 +60,15 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non Ak_tmp = Ak.conj() # Rotate the Ak wavefunction using the representation matrices - rot_Ak[:, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp) @ (Dvv.transpose(0, 2, 1)[None, ...]) - + if tda : + rot_Ak[:, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp) @ (Dvv.transpose(0, 2, 1)[None, ...]) + else : + ## rotate the resonant part + rot_Ak[:, 0, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp[:,0,...]) @ (Dvv.transpose(0, 2, 1)[None, ...]) + # Rotate the anti-resonant part + Dvv = dmats[:, :nv, :nv] + Dcc = dmats[idx_k_minus_q, nv:, nv:].conj() + rot_Ak[:, 1, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp[:,1,...]) @ (Dvv.transpose(0, 2, 1)[None, ...]) + # done return the rotate wfc return rot_Ak diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index 9ceb9fa8..07359cfc 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -90,6 +90,7 @@ def from_db_file(cls,lattice,filename='ndb.BS_diago_Q1',folder='.',Load_WF=True, Set `Read_WF=False` to avoid reading eigenvectors for faster IO and memory efficiency. If neigs < 0 ; all eigen values (vectors) are loaded or else first neigs are loaded + " In case of non-TDA, we load right eigenvectors. """ path_filename = os.path.join(folder,filename) if not os.path.isfile(path_filename): @@ -268,14 +269,14 @@ def get_Akcv(self): sort_idx = bs_table0*nc*nv + bs_table2*nv + bs_table1 # eig_wfcs_returned[:,sort_idx] = eig_wfcs[...,:table_len] - eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) - # # check if this is coupling . - # if eig_wfcs.shape[-1]//table_len == 2: - # eig_wfcs_returned[:,sort_idx+table_len] = eig_wfcs[...,table_len:] - # eig_wfcs_returned = eig_wfcs_returned.reshape(-1,2,nk,nc,nv) - # else : - # eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) + if eig_wfcs.shape[-1]//table_len == 2: + eig_wfcs_returned[:,sort_idx+table_len] = eig_wfcs[...,table_len:] + # NM : Note that here v and c are inverted i.e + # psi_S = Akcv * phi_v(r_e) * phi_c^*(r_h) + eig_wfcs_returned = eig_wfcs_returned.reshape(-1,2,nk,nc,nv) + else : + eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) # self.Akcv = eig_wfcs_returned return self.Akcv From 5289dab61dce1f7b800d3dee978ca833c88669a2 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 25 Mar 2025 17:37:39 +0100 Subject: [PATCH 32/73] remove assert to allow non TDA --- yambopy/dbs/excitondb.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index 07359cfc..a5e12482 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -255,8 +255,6 @@ def get_Akcv(self): # Make sure nc * nv * nk = BS_TABLE length table_len = nk*nv*nc assert table_len == self.table.shape[0], "BS_TABLE length not equal to nc * nv * nk" - - assert eig_wfcs.shape[-1]//table_len == 1, "rearranged_Akcv works only in TDA" # v_min = np.min(self.table[:,1]) c_min = np.min(self.table[:,2]) From c58b835350be215c2fba517d0196c53013fb50c7 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 1 Apr 2025 10:04:42 +0200 Subject: [PATCH 33/73] fix a tiny bug --- yambopy/bse/exciton_spin.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index a5428556..c508e768 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -169,7 +169,7 @@ def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, def get_spinvals(spin_matrix, eigenvalues, atol=1e-3, rtol=1e-3): - degen_idx = find_degeneracy_evs(eigenvalues) + degen_idx = find_degeneracy_evs(eigenvalues,atol=atol, rtol=rtol) spins = [] for id in degen_idx: w = np.linalg.eigvals(spin_matrix[id,:][:,id]) From 4a12097ca937840a7609df3214c5726b4cc77864 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 1 Apr 2025 10:28:21 +0200 Subject: [PATCH 34/73] slightly lower the defauly degen threshould --- yambopy/bse/exciton_spin.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index c508e768..8f40188d 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -72,7 +72,7 @@ def compute_exciton_spin(lattice, excdb, wfdb, elec_sz, contribution='b',diagona def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, - nstates=-1, contribution='b', degen_tol = 1e-3, + nstates=-1, contribution='b', degen_tol = 1e-2, sz=0.5 * np.array([[1, 0], [0, -1]]), return_dbs_and_spin=True): """ @@ -99,7 +99,7 @@ def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, - `'e'`: Electron spin only. - `'h'`: Hole spin only. degen_tol : float, optional - Degeneracy tolerance for excitons in eV. Default is `1e-3` eV. + Degeneracy tolerance for excitons in eV. Default is `1e-2` eV. sz : ndarray, optional Spin-z operator matrix in the basis of spinor wavefunctions. Default is `0.5 * np.array([[1, 0], [0, -1]])`. From 9298813029d60a02578e0215986ea899901c12b0 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 1 Apr 2025 22:43:19 +0200 Subject: [PATCH 35/73] bug in cube file format --- yambopy/io/cubetools.py | 1 + 1 file changed, 1 insertion(+) diff --git a/yambopy/io/cubetools.py b/yambopy/io/cubetools.py index d5f11100..5b7d0f2d 100644 --- a/yambopy/io/cubetools.py +++ b/yambopy/io/cubetools.py @@ -35,4 +35,5 @@ def write_cube(filename, data, lat_vec, atom_pos, atomic_num, origin=np.zeros(3) for iz in range(data.shape[2]): cube.write("%.6f " %(data[ix,iy,iz])) if iz % 6 == 5 : cube.write("\n") + cube.write("\n") From 158c96938372680afb77ec22278d69cc96ee3bd5 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 1 Apr 2025 22:59:22 +0200 Subject: [PATCH 36/73] fix bug in cubetools --- yambopy/io/cubetools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yambopy/io/cubetools.py b/yambopy/io/cubetools.py index 5b7d0f2d..934619e2 100644 --- a/yambopy/io/cubetools.py +++ b/yambopy/io/cubetools.py @@ -11,10 +11,10 @@ def write_cube(filename, data, lat_vec, atom_pos, atomic_num, origin=np.zeros(3) atomic_num: atomic numbers of elements """ ## atomic positions must be in [0,1) - apos_crys = atom_pos@lat_vec.T + apos_crys = atom_pos@np.linalg.inv(lat_vec.T) apos_crys = apos_crys-np.floor(apos_crys) apos_crys = (apos_crys+1e-6)%1 - apos_new = apos_crys@np.linalg.inv(lat_vec.T) + apos_new = apos_crys@lat_vec.T with open(filename, "w") as cube: cube.write("# Cube file generated using YamboPy\n") From 664f2481c0fea82962ba9e91acb62aed6e5c7c9c Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 8 Apr 2025 16:02:29 +0200 Subject: [PATCH 37/73] First version for real space ex_wf with phase information. --- yambopy/bse/realSpace_excitonwf.py | 274 +++++++++++++++++++++++++++++ yambopy/dbs/wfdb.py | 2 +- 2 files changed, 275 insertions(+), 1 deletion(-) create mode 100644 yambopy/bse/realSpace_excitonwf.py diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py new file mode 100644 index 00000000..133b9d54 --- /dev/null +++ b/yambopy/bse/realSpace_excitonwf.py @@ -0,0 +1,274 @@ +### Compute real space exction wavefunction when hole/electron is fixed. +import numpy as np +from scipy.spatial import KDTree +from yambopy.kpoints import build_ktree, find_kpt +from yambopy.dbs.wfdb import su2_mat +from yambopy.dbs.excitondb import YamboExcitonDB +from yambopy.dbs.latticedb import YamboLatticeDB +from yambopy.dbs.wfdb import YamboWFDB +from yambopy.io.cubetools import write_cube +import os + +def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, + fix_particle='h', supercell=[1,1,1], wfcCutoffRy=-1, block_size=256): + # NM : please note that this only computes Akcv * psi_{kc}(r_e) * (psi_{k-Q,v}(r_h))^* + # For density, one must compute the abosulte value. + # + # + # Akcv, (Nstates,k,c,v) : Exciton wfc. No checking is done internally about the shape. + # Qpt : Qpt point of exciton in crystal coordinates + # wfcdb : wf db + # band indices used in bse [nb1, nb2]. + # nb1 and nb2 are same indices used in yambo input i.e % BSEBands nb1 | nb2 % in yambo input + # fixed_postion : postion of electron or hole in crystal coordinates + # fix_particle : 'e' : electron positon is fixed and hole density is computed + # 'h' : hole is fixed and the electronic dentisy is computed (default) + # supercell : size of supercell list or array of 3 intergers. If even integergs + # are given, we add +1 and turn it to odd. + # wfcCutoffRy : Cutoff for Wfc in Rydberg units, default = -1 (full wfc is cutoff is used) + # block_size : is a postive integer, the default is 256 which is very good but uses more memory. + # decrease it when you run into memory issues + + # Convert them to + for i in range(3): + if supercell[i]%2 == 0: + print('Warning : Even supercell given, so increasing' + + ' supercell size along %d direction by 1'%(i+1)) + supercell[i] = supercell[i] +1 + # + # + fix_particle = fix_particle.lower() + bse_bnds = [min(bse_bnds)-1,max(bse_bnds)] + assert bse_bnds[0] >= wfcdb.min_bnd, \ + "%d is used in bse but not found in wfcdb, load more wfcs" %(bse_bnds[0]+1) + assert bse_bnds[1] <= wfcdb.min_bnd + wfcdb.nbands, \ + "%d is used in bse but not found in wfcdb, load more wfcs" %(wfcdb.min_bnd + wfcdb.nbands) + bse_bnds = np.array(bse_bnds,dtype=int)-wfcdb.min_bnd + + nstates, nk, nc, nv = Akcv.shape + + kpt_idx = wfcdb.ydb.kpoints_indexes + sym_idx = wfcdb.ydb.symmetry_indexes + nkBZ = len(sym_idx) + assert nc + nv == bse_bnds[1]-bse_bnds[0], "Band mismatch" + assert nk == nkBZ, "kpoint mismatch" + # + + ## Bring the positon of hole/electron in the centre unit cell, i.e make in between [0,1) + fixed_postion = np.array(fixed_postion) - np.floor(fixed_postion) + fixed_postion = (fixed_postion + 1e-6)%1 + fixed_postion = fixed_postion - 1e-6 + + fixed_postion += np.array(supercell)//2 + + + hole_bnds = [bse_bnds[0],bse_bnds[0]+nv] + elec_bnds = [bse_bnds[0]+nv,bse_bnds[1]] + + lat_vec = wfcdb.ydb.lat.T + blat = np.linalg.inv(lat_vec) + gvecs_iBZ_idx = [] + fft_box = np.zeros(3,dtype=int) + # + + for ik in range(len(wfcdb.gvecs)): + idx_gvecs_tmp = np.arange(wfcdb.ngvecs[ik],dtype=int) + if wfcCutoffRy > 0: + tmp_gvecs = np.linalg.norm((wfcdb.gvecs[ik, :wfcdb.ngvecs[ik], :] + + wfcdb.kpts_iBZ[ik][None,:])@blat,axis=-1) + idx_tmp = tmp_gvecs < np.sqrt(wfcCutoffRy) + idx_gvecs_tmp = idx_gvecs_tmp[idx_tmp].copy() + # + gvecs_iBZ_idx.append(idx_gvecs_tmp) + ## Get the fft box + min_fft_idx = np.min(wfcdb.gvecs[ik, :wfcdb.ngvecs[ik], :][idx_gvecs_tmp] , axis=0) + max_fft_idx = np.max(wfcdb.gvecs[ik, :wfcdb.ngvecs[ik], :][idx_gvecs_tmp] , axis=0) + assert np.min(max_fft_idx) >= 0 and np.max(min_fft_idx) < 0, "Invalid G-vectors" + for i in range(3): + fft_box[i] = max([fft_box[i], max_fft_idx[i] - min_fft_idx[i] + 3]) + + # Compute nstates, nk, Nx, Ny, Nz object + print("FFT Box : ",fft_box[0], fft_box[1], fft_box[2]) + # + ktree = build_ktree(wfcdb.kBZ) + # + nspinorr = wfcdb.nspinor + exe_wfc_real = np.zeros((nstates, nspinorr, np.prod(supercell), + fft_box[0], fft_box[1], fft_box[2]),dtype=np.complex64) + Lx = np.arange(supercell[0],dtype=int) + Ly = np.arange(supercell[1],dtype=int) + Lz = np.arange(supercell[2],dtype=int) + Lsupercells = np.zeros((supercell[0],supercell[1],supercell[2],3),dtype=int) + Lsupercells[...,0], Lsupercells[...,1], Lsupercells[...,2] = np.meshgrid(Lx,Ly,Lz,indexing='ij') + # + FFTxx = np.fft.fftfreq(fft_box[0]) + FFTyy = np.fft.fftfreq(fft_box[1]) + FFTzz = np.fft.fftfreq(fft_box[2]) + FFTxx = FFTxx - np.floor(FFTxx) + FFTyy = FFTyy - np.floor(FFTyy) + FFTzz = FFTzz - np.floor(FFTzz) + FFFboxs = np.zeros((fft_box[0],fft_box[1],fft_box[2],3)) + FFFboxs[...,0], FFFboxs[...,1], FFFboxs[...,2] = np.meshgrid(FFTxx,FFTyy,FFTzz,indexing='ij') + # + exe_tmp_wf = np.zeros((nstates, nspinorr, min(nk,block_size) , + fft_box[0], fft_box[1], fft_box[2]),dtype=np.complex64) + # + exp_tmp_kL = np.zeros((min(nk,block_size) ,supercell[0], + supercell[1],supercell[2]),dtype=np.complex64) + + nblks = nk//block_size + nrem = nk%block_size + if nrem > 0: nblks = nblks+1 + # + for ibk in range(nblks): + ikstart = ibk*block_size + ikstop = min(ikstart + block_size,nk) + for ik in range(ikstart,ikstop): + ## First get the electronic wfcs + ik_ibz = kpt_idx[ik] + isym = sym_idx[ik] + wfc_tmp, gvecs_tmp = wfcdb.get_iBZ_wf(ik_ibz) + wfc_tmp = wfc_tmp[:,elec_bnds[0]:elec_bnds[1],:,gvecs_iBZ_idx[ik_ibz]] + gvecs_tmp = gvecs_tmp[gvecs_iBZ_idx[ik_ibz]] + kvec = wfcdb.get_iBZ_kpt(ik_ibz) + # get the rotated wf + if isym != 0: + sym_mat = wfcdb.ydb.sym_car[isym] + time_rev = (isym >= len(wfcdb.ydb.sym_car + ) / (1 + int(np.rint(wfcdb.ydb.time_rev)))) + kvec, wfc_tmp, gvecs_tmp = wfcdb.apply_symm( + kvec, wfc_tmp, gvecs_tmp, time_rev, sym_mat) + + kelec = kvec + wfc_elec = wfc_tmp + gvecs_elec = gvecs_tmp + + ## Do the same and get hole wfc + ikhole = find_kpt(ktree, kelec-Qpt) + ik_ibz = kpt_idx[ikhole] + isym = sym_idx[ikhole] + wfc_tmp, gvecs_tmp = wfcdb.get_iBZ_wf(ik_ibz) + wfc_tmp = wfc_tmp[:,hole_bnds[0]:hole_bnds[1],:,gvecs_iBZ_idx[ik_ibz]] + gvecs_tmp = gvecs_tmp[gvecs_iBZ_idx[ik_ibz]] + kvec = wfcdb.get_iBZ_kpt(ik_ibz) + # get the rotated wf + if isym != 0: + sym_mat = wfcdb.ydb.sym_car[isym] + time_rev = (isym >= len(wfcdb.ydb.sym_car) / (1 + int(np.rint(wfcdb.ydb.time_rev)))) + kvec, wfc_tmp, gvecs_tmp = wfcdb.apply_symm(kvec, wfc_tmp, gvecs_tmp, time_rev, sym_mat) + # + khole = -kvec + wfc_hole = wfc_tmp.conj() + gvecs_hole = -gvecs_tmp + + if fix_particle == 'h': + fx_kvec = khole + fx_wfc = wfc_hole + fx_gvec = gvecs_hole + # + ft_kvec = kelec + ft_wfc = wfc_elec + ft_gvec = gvecs_elec + else : + ft_kvec = khole + ft_wfc = wfc_hole + ft_gvec = gvecs_hole + # + fx_kvec = kelec + fx_wfc = wfc_elec + fx_gvec = gvecs_elec + # compute + ## Fix compute the fixed particle wfc in real space. + ## NM : Donot perform FFT as we only need it for one point. + exp_fx = np.exp(2*np.pi*1j*((fx_gvec + fx_kvec[None,:])@fixed_postion)) + fx_wfc *= exp_fx[None,None,None,:] + fx_wfc = np.sum(fx_wfc,axis=-1) #(spin,bnd,spinor) + ## for now only nspin = 1 works. + ns1, nbndc, nspinorr, ng = ft_wfc.shape + #if ft_ikpt not in prev_ikpts: + ft_wfcr = wfcdb.to_real_space(ft_wfc.reshape(-1,nspinorr,ng),ft_gvec, grid=fft_box) + ft_wfcr = ft_wfcr.reshape(ns1,nbndc,nspinorr,fft_box[0],fft_box[1],fft_box[2]) + exp_kx_r = np.exp(2*np.pi*1j*FFFboxs.reshape(-1,3)@ft_kvec) + # + if fix_particle == 'h': + np.einsum('ncv,vx,cxijk->nxijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], + optimize=True,out=exe_tmp_wf[:,:,ik-ikstart]) + else : + np.einsum('ncv,cx,vxijk->nxijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], + optimize=True,out=exe_tmp_wf[:,:,ik-ikstart]) + # + exe_tmp_wf[:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3]) + exp_tmp_kL[ik-ikstart] = np.exp(1j*2*np.pi*np.einsum('...x,x->...',Lsupercells,ft_kvec)) + ## perform gemm operation + exe_wfc_real[...] += (exp_tmp_kL.reshape(len(exp_tmp_kL),-1)[:(ikstop-ikstart)].T[None,...] @ + exe_tmp_wf.reshape(nstates*nspinorr,-1,np.prod(fft_box))[:,:(ikstop-ikstart)]).reshape( + nstates,nspinorr,np.prod(supercell),fft_box[0], fft_box[1], fft_box[2]) + + exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr,supercell[0],supercell[1],supercell[2], + fft_box[0], fft_box[1], fft_box[2]) + # + exe_wfc_real = exe_wfc_real.transpose(0,1,2,5,3,6,4,7).reshape(nstates,nspinorr, + supercell[0]*fft_box[0], supercell[1]*fft_box[1], + supercell[2]*fft_box[2])/np.prod(supercell) + + # compute postioon of fixed particle in cart units + fixed_postion_cc = lat_vec@fixed_postion + Lsupercells = Lsupercells.reshape(-1,3)#/np.array(supercell)[None,:] + Lsupercells = Lsupercells@lat_vec.T + atom_pos = Lsupercells[:,None,:] + wfcdb.ydb.car_atomic_positions[None,:,:] + atom_pos = np.append(atom_pos.reshape(-1,3),fixed_postion_cc[None,:],axis=0) + supercell_latvecs = lat_vec*np.array(supercell)[None,:] + ## Make atomic numbers + atom_nums = np.zeros(len(atom_pos),dtype=wfcdb.ydb.atomic_numbers.dtype) + attmp = atom_nums[:-1].reshape(-1,len(wfcdb.ydb.atomic_numbers)) + attmp[...]= wfcdb.ydb.atomic_numbers[None,:] + atom_nums[-1] = 200 + # + return supercell_latvecs,atom_nums,atom_pos,exe_wfc_real + +def compute_exc_wfc_real(path='.', bse_dir='SAVE', iqpt=1, nstates=[1], + fixed_postion=[0,0,0], fix_particle='h', aveg=True, supercell=[1,1,1], + wfcCutoffRy=-1, phase=False, block_size=256): + # + lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) + filename = 'ndb.BS_diago_Q%d' % (iqpt) + excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, + folder=os.path.join(path, bse_dir), + Load_WF=True, neigs=max(nstates)) + # Load the wavefunction database + wfdb = YamboWFDB(path=path, latdb=lattice, + bands_range=[np.min(excdb.table[:, 1]) - 1, + np.max(excdb.table[:, 2])]) + # + Akcv = excdb.get_Akcv()[min(nstates)-1:max(nstates)] + excQpt = excdb.car_qpoint + # + # Convert the q-point to crystal coordinates + Qpt = wfdb.ydb.lat @ excQpt + + sc_latvecs, atom_nums, atom_pos, real_wfc = ex_wf2Real(Akcv, Qpt, wfdb, [np.min(excdb.table[:, 1]), + np.max(excdb.table[:, 2])], fixed_postion=fixed_postion, + fix_particle=fix_particle, supercell=supercell, + wfcCutoffRy=wfcCutoffRy, block_size=block_size) + # + nstates_range = np.arange(nstates[0],nstates[1],dtype=int) + density = np.abs(real_wfc)**2 + if phase: + phase = np.sign(real_wfc.real) #np.sign(np.angle(real_wfc)) + density *= phase + if aveg: + real_wfc = np.sum(density,axis=(0,1)) + real_wfc = real_wfc/np.max(np.abs(real_wfc)) + write_cube('exe_wf_avg_%d-%d.cube' %(nstates[0],nstates[1]), + real_wfc, sc_latvecs, atom_pos, atom_nums, header='Real space exciton wavefunction') + else: + real_wfc = np.sum(density,axis=1) + for i in range(len(real_wfc)): + real_wfc1 = real_wfc[i]/np.max(np.abs(real_wfc[i])) + write_cube('exe_wf_%d.cube' %nstates[i], real_wfc1, sc_latvecs, + atom_pos, atom_nums, header='Real space exciton wavefunction') + + + +#if __name__ == "__main__": + diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 239debdc..4b55c723 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -442,7 +442,7 @@ def to_real_space(self, wfc_tmp, gvec_tmp, grid=[]): grid = self.fft_box cel_vol = abs(np.linalg.det(self.ydb.lat.T)) - tmp_wfc = np.zeros((self.nspin, self.nspinor, grid[0], grid[1], grid[2]), dtype=wfc_tmp.dtype) + tmp_wfc = np.zeros((len(wfc_tmp), self.nspinor, grid[0], grid[1], grid[2]), dtype=wfc_tmp.dtype) Nx_vals = np.where(gvec_tmp[:, 0] >= 0, gvec_tmp[:, 0], gvec_tmp[:, 0] + grid[0]) Ny_vals = np.where(gvec_tmp[:, 1] >= 0, gvec_tmp[:, 1], gvec_tmp[:, 1] + grid[1]) From 858353d210ad1d9ea0d3267698de3d2f8414eeea Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 8 Apr 2025 18:50:19 +0200 Subject: [PATCH 38/73] Fix bug in case of spinors in ex_wf --- yambopy/bse/realSpace_excitonwf.py | 37 ++++++++++++++++++------------ 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 133b9d54..cfcfb9fb 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -93,7 +93,7 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, ktree = build_ktree(wfcdb.kBZ) # nspinorr = wfcdb.nspinor - exe_wfc_real = np.zeros((nstates, nspinorr, np.prod(supercell), + exe_wfc_real = np.zeros((nstates, nspinorr, nspinorr, np.prod(supercell), fft_box[0], fft_box[1], fft_box[2]),dtype=np.complex64) Lx = np.arange(supercell[0],dtype=int) Ly = np.arange(supercell[1],dtype=int) @@ -110,7 +110,7 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, FFFboxs = np.zeros((fft_box[0],fft_box[1],fft_box[2],3)) FFFboxs[...,0], FFFboxs[...,1], FFFboxs[...,2] = np.meshgrid(FFTxx,FFTyy,FFTzz,indexing='ij') # - exe_tmp_wf = np.zeros((nstates, nspinorr, min(nk,block_size) , + exe_tmp_wf = np.zeros((nstates, nspinorr, nspinorr, min(nk,block_size) , fft_box[0], fft_box[1], fft_box[2]),dtype=np.complex64) # exp_tmp_kL = np.zeros((min(nk,block_size) ,supercell[0], @@ -191,23 +191,23 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, exp_kx_r = np.exp(2*np.pi*1j*FFFboxs.reshape(-1,3)@ft_kvec) # if fix_particle == 'h': - np.einsum('ncv,vx,cxijk->nxijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], - optimize=True,out=exe_tmp_wf[:,:,ik-ikstart]) + np.einsum('ncv,vy,cxijk->nxyijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], + optimize=True,out=exe_tmp_wf[:,:,:,ik-ikstart]) else : - np.einsum('ncv,cx,vxijk->nxijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], - optimize=True,out=exe_tmp_wf[:,:,ik-ikstart]) + np.einsum('ncv,cx,vyijk->nxyijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], + optimize=True,out=exe_tmp_wf[:,:,:,ik-ikstart]) # - exe_tmp_wf[:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3]) + exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] exp_tmp_kL[ik-ikstart] = np.exp(1j*2*np.pi*np.einsum('...x,x->...',Lsupercells,ft_kvec)) ## perform gemm operation exe_wfc_real[...] += (exp_tmp_kL.reshape(len(exp_tmp_kL),-1)[:(ikstop-ikstart)].T[None,...] @ - exe_tmp_wf.reshape(nstates*nspinorr,-1,np.prod(fft_box))[:,:(ikstop-ikstart)]).reshape( - nstates,nspinorr,np.prod(supercell),fft_box[0], fft_box[1], fft_box[2]) + exe_tmp_wf.reshape(nstates*nspinorr**2,-1,np.prod(fft_box))[:,:(ikstop-ikstart)]).reshape( + nstates,nspinorr,nspinorr,np.prod(supercell),fft_box[0], fft_box[1], fft_box[2]) - exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr,supercell[0],supercell[1],supercell[2], + exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr**2,supercell[0],supercell[1],supercell[2], fft_box[0], fft_box[1], fft_box[2]) # - exe_wfc_real = exe_wfc_real.transpose(0,1,2,5,3,6,4,7).reshape(nstates,nspinorr, + exe_wfc_real = exe_wfc_real.transpose(0,1,2,5,3,6,4,7).reshape(nstates,nspinorr,nspinorr, supercell[0]*fft_box[0], supercell[1]*fft_box[1], supercell[2]*fft_box[2])/np.prod(supercell) @@ -251,21 +251,28 @@ def compute_exc_wfc_real(path='.', bse_dir='SAVE', iqpt=1, nstates=[1], fix_particle=fix_particle, supercell=supercell, wfcCutoffRy=wfcCutoffRy, block_size=block_size) # + # nstates_range = np.arange(nstates[0],nstates[1],dtype=int) density = np.abs(real_wfc)**2 + + if fix_particle == 'h': name_file = 'electron' + else: name_file = 'hole' + if real_wfc.shape[1] != 1: + print("phase plot only works for nspin = 1 and nspinor == 1") + phase = False if phase: phase = np.sign(real_wfc.real) #np.sign(np.angle(real_wfc)) density *= phase if aveg: - real_wfc = np.sum(density,axis=(0,1)) + real_wfc = np.sum(density,axis=(0,1,2)) real_wfc = real_wfc/np.max(np.abs(real_wfc)) - write_cube('exe_wf_avg_%d-%d.cube' %(nstates[0],nstates[1]), + write_cube('exe_wf_avg_%s_%d-%d.cube' %(name_file,nstates[0],nstates[1]), real_wfc, sc_latvecs, atom_pos, atom_nums, header='Real space exciton wavefunction') else: - real_wfc = np.sum(density,axis=1) + real_wfc = np.sum(density,axis=(1,2)) for i in range(len(real_wfc)): real_wfc1 = real_wfc[i]/np.max(np.abs(real_wfc[i])) - write_cube('exe_wf_%d.cube' %nstates[i], real_wfc1, sc_latvecs, + write_cube('exe_wf_%s_%d.cube' %(name_file,nstates[i]), real_wfc1, sc_latvecs, atom_pos, atom_nums, header='Real space exciton wavefunction') From a79af2769ce7e6ec575ad23819c2768ba275c8f4 Mon Sep 17 00:00:00 2001 From: Murali Date: Wed, 9 Apr 2025 10:51:24 +0200 Subject: [PATCH 39/73] improve memory usage --- yambopy/bse/realSpace_excitonwf.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index cfcfb9fb..a5f7e954 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -200,16 +200,25 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] exp_tmp_kL[ik-ikstart] = np.exp(1j*2*np.pi*np.einsum('...x,x->...',Lsupercells,ft_kvec)) ## perform gemm operation - exe_wfc_real[...] += (exp_tmp_kL.reshape(len(exp_tmp_kL),-1)[:(ikstop-ikstart)].T[None,...] @ - exe_tmp_wf.reshape(nstates*nspinorr**2,-1,np.prod(fft_box))[:,:(ikstop-ikstart)]).reshape( - nstates,nspinorr,nspinorr,np.prod(supercell),fft_box[0], fft_box[1], fft_box[2]) + total_gemms_t = nstates*nspinorr**2 + exp_tmp_kL_tmp = exp_tmp_kL.reshape(len(exp_tmp_kL),-1)[:(ikstop-ikstart)].T + exe_tmp_wf_tmp = exe_tmp_wf.reshape(nstates,nspinorr,nspinorr,-1,np.prod(fft_box)) + exe_tmp_wf_tmp = exe_tmp_wf_tmp[...,:(ikstop-ikstart),:] + # + for igemms in range(total_gemms_t): + ii, jj, kk = np.unravel_index(igemms, (nstates,nspinorr,nspinorr)) + # NM : It is not nice to create an large temporary array again. but numpy does support + # C += A@B call like that blas has. + exe_wfc_real[ii, jj, kk ] += (exp_tmp_kL_tmp @ exe_tmp_wf_tmp[ii, jj, kk ]).reshape( + np.prod(supercell),fft_box[0], fft_box[1], fft_box[2]) exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr**2,supercell[0],supercell[1],supercell[2], fft_box[0], fft_box[1], fft_box[2]) # exe_wfc_real = exe_wfc_real.transpose(0,1,2,5,3,6,4,7).reshape(nstates,nspinorr,nspinorr, - supercell[0]*fft_box[0], supercell[1]*fft_box[1], - supercell[2]*fft_box[2])/np.prod(supercell) + supercell[0]*fft_box[0], supercell[1]*fft_box[1], + supercell[2]*fft_box[2]) + exe_wfc_real *= (1.0/np.prod(supercell)) # compute postioon of fixed particle in cart units fixed_postion_cc = lat_vec@fixed_postion From 2734892021e660814938cc40c840014154b6f1d1 Mon Sep 17 00:00:00 2001 From: Murali Date: Wed, 9 Apr 2025 11:02:48 +0200 Subject: [PATCH 40/73] add a small check --- yambopy/bse/realSpace_excitonwf.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index a5f7e954..052178c1 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -29,6 +29,10 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # block_size : is a postive integer, the default is 256 which is very good but uses more memory. # decrease it when you run into memory issues + if block_size < 1: + print('Warning: Wrong block_size. setting to 1') + block_size = 1 + # # Convert them to for i in range(3): if supercell[i]%2 == 0: From 7d0cd3e9d6e39588a842b4a4961d9d323c8d5452 Mon Sep 17 00:00:00 2001 From: Murali Date: Wed, 9 Apr 2025 23:01:55 +0200 Subject: [PATCH 41/73] remove unnecessary imports --- yambopy/bse/realSpace_excitonwf.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 052178c1..e7940369 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -1,8 +1,6 @@ ### Compute real space exction wavefunction when hole/electron is fixed. import numpy as np -from scipy.spatial import KDTree from yambopy.kpoints import build_ktree, find_kpt -from yambopy.dbs.wfdb import su2_mat from yambopy.dbs.excitondb import YamboExcitonDB from yambopy.dbs.latticedb import YamboLatticeDB from yambopy.dbs.wfdb import YamboWFDB From c52a3e4a6354d884b2a875dc0d070c82d70b60a8 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 10 Apr 2025 11:14:28 +0200 Subject: [PATCH 42/73] add comments --- yambopy/bse/realSpace_excitonwf.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index e7940369..8d7d64a1 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -27,6 +27,14 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # block_size : is a postive integer, the default is 256 which is very good but uses more memory. # decrease it when you run into memory issues + ## ! Natoms_in_supercell = (natom_in_unit_call * Nsupercell) + 1 (+1 due to hole/electron) + # Outputs : + # supercell_latvecs : New lattice vectors for the big supercell (3,3) + # atom_nums : atomic numbers of atoms in supercell (Natoms_in_supercell) + # atom_pos : atomic positons in cart units for atoms in supercell. (Natoms_in_supercell,3) + # exe_wfc_real: complex array (nstates,nspinor_electron,nspinor_hole, FFTx, FFTy, FFTz). + # Please note that you must take absoulte square and contract the spinor dimension to the density. + # if block_size < 1: print('Warning: Wrong block_size. setting to 1') block_size = 1 From e9e51fed9f1bf1bec18ceb58add3c540197f2a2d Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 13 Apr 2025 23:20:34 +0200 Subject: [PATCH 43/73] finalize real space exciton wfc --- yambopy/bse/realSpace_excitonwf.py | 310 ++++++++++++++++++++--------- yambopy/dbs/excitondb.py | 75 ++++++- 2 files changed, 290 insertions(+), 95 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 8d7d64a1..31d9d202 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -4,36 +4,145 @@ from yambopy.dbs.excitondb import YamboExcitonDB from yambopy.dbs.latticedb import YamboLatticeDB from yambopy.dbs.wfdb import YamboWFDB -from yambopy.io.cubetools import write_cube +from tqdm import tqdm import os -def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, - fix_particle='h', supercell=[1,1,1], wfcCutoffRy=-1, block_size=256): - # NM : please note that this only computes Akcv * psi_{kc}(r_e) * (psi_{k-Q,v}(r_h))^* - # For density, one must compute the abosulte value. - # - # - # Akcv, (Nstates,k,c,v) : Exciton wfc. No checking is done internally about the shape. - # Qpt : Qpt point of exciton in crystal coordinates - # wfcdb : wf db - # band indices used in bse [nb1, nb2]. - # nb1 and nb2 are same indices used in yambo input i.e % BSEBands nb1 | nb2 % in yambo input - # fixed_postion : postion of electron or hole in crystal coordinates - # fix_particle : 'e' : electron positon is fixed and hole density is computed - # 'h' : hole is fixed and the electronic dentisy is computed (default) - # supercell : size of supercell list or array of 3 intergers. If even integergs - # are given, we add +1 and turn it to odd. - # wfcCutoffRy : Cutoff for Wfc in Rydberg units, default = -1 (full wfc is cutoff is used) - # block_size : is a postive integer, the default is 256 which is very good but uses more memory. +## Usage +""" +import numpy as np +from yambopy.dbs.excitondb import YamboExcitonDB +from yambopy.dbs.latticedb import YamboLatticeDB +from yambopy.dbs.wfdb import YamboWFDB + +iqpt = 1 # qpt index of exciton + +# load lattice db +lattice = YamboLatticeDB.from_db_file(os.path.join('.', '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(path, bse_dir) + neigs = 20) + +#Load the wavefunction database +wfdb = YamboWFDB(path=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 +# 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, fixed_postion=[0,0,0], supercell=[1,1,1], + degen_tol=0.01, wfcCutoffRy=80, fix_particle='h') + +## .cube will be dumped and use vesta to visualize it ! +""" + +## + +def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, + fix_particle='h', supercell=[1,1,1], wfcCutoffRy=-1, + block_size=256): + """ + Compute real-space exciton wavefunction when hole/electron is fixed. + + This is the main interface function that handles both resonant and anti-resonant parts + of the exciton wavefunction in the case of non-TDA calculations. + + Args: + Akcv (numpy.ndarray): Exciton wavefunction coefficients with shape: + - (Nstates,k,c,v) for TDA + - (Nstates,2,k,c,v) for non-TDA (2 for resonant/anti-resonant) + Qpt (numpy.ndarray): Q-point of exciton in crystal coordinates + wfcdb (YamboWFDB): Wavefunction database + bse_bnds (list): Band range used in BSE [min_band, max_band] (python indexing) + fixed_postion (list): Position of fixed particle in crystal coordinates + fix_particle (str): 'e' to fix electron, 'h' to fix hole (default) + supercell (list): Supercell dimensions [nx,ny,nz] + wfcCutoffRy (float): Wavefunction cutoff in Rydberg (-1 for no cutoff) + block_size (int): Block size for memory-efficient computation. + ## choosing lowe block_size will slight lower the memory requirment but also less faster + + Returns: + tuple: (supercell_latvecs, atom_nums, atom_pos, exe_wfc_real) + - supercell_latvecs: Supercell lattice vectors (3,3) + - atom_nums: Atomic numbers. + - atom_pos: Atomic positions in cartisian units + - exe_wfc_real: Real-space exciton wavefunction (nstates, nspinor_electron, nspinor_hole, + Nx_grid, Ny_grid, Nz_grid) + """ + if len(Akcv.shape) == 5: + #nonTDA + ## first the resonat part + supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ + ex_wf2Real_kernel(Akcv[:,0], Qpt, wfcdb, bse_bnds, fixed_postion, + fix_particle=fix_particle, supercell=supercell, + wfcCutoffRy=wfcCutoffRy, block_size=block_size, + ares=False, out_res=None) + # add to antiresonant part + supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ + ex_wf2Real_kernel(Akcv[:,1], Qpt, wfcdb, bse_bnds, fixed_postion, + fix_particle=fix_particle, supercell=supercell, + wfcCutoffRy=wfcCutoffRy, block_size=block_size, + ares=True, out_res=exe_wfc_real) + else: + assert len(Akcv.shape) == 4, "wrong Akcv dimensions" + supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ + ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, + fix_particle=fix_particle, supercell=supercell, + wfcCutoffRy=wfcCutoffRy, block_size=block_size, + ares=False, out_res=None) + + return supercell_latvecs,atom_nums,atom_pos,exe_wfc_real + + + +def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, + fix_particle='h', supercell=[1,1,1], wfcCutoffRy=-1, + block_size=256, ares=False, out_res=None): + """ + Core kernel function for computing real-space exciton wavefunction. + + Computes either: + - Akcv * psi_{kc}(r_e) * (psi_{k-Q,v}(r_h))^* (resonant part) + - Akcv * psi_{kv}(r_e) * (psi_{k-Q,c}(r_h))^* (anti-resonant part) + + Note: For density, one must compute the absolute value squared. + + Args: + Akcv (numpy.ndarray): Exciton coefficients [nstates,nk,nc,nv] + Qpt (numpy.ndarray): Q-point in crystal coordinates [3] + wfcdb (YamboWFDB): Wavefunction database + bse_bnds (list): BSE band range used in bse [nb1, nb2]. fortran indexing. i.e index starts from 1. + i.e nb1 and nb2 are same indices used in yambo input i.e + % BSEBands nb1 | nb2 % in yambo input + fixed_postion (list): Fixed particle position in crystal coords [3] + fix_particle (str): 'e'=fix electron, 'h'=fix hole (default) + supercell (list): Supercell dimensions [nx,ny,nz] + wfcCutoffRy (float): Wavefunction cutoff in Rydberg (-1= full cutoff) + block_size (int): Memory block size for computation. is a postive integer, + the default is 256 which is very good but uses more memory. # decrease it when you run into memory issues + ares (bool): If True, compute anti-resonant part + out_res (numpy.ndarray): If provided and ares=True, adds to this array + and is returned instead of internally creating. - ## ! Natoms_in_supercell = (natom_in_unit_call * Nsupercell) + 1 (+1 due to hole/electron) - # Outputs : - # supercell_latvecs : New lattice vectors for the big supercell (3,3) - # atom_nums : atomic numbers of atoms in supercell (Natoms_in_supercell) - # atom_pos : atomic positons in cart units for atoms in supercell. (Natoms_in_supercell,3) - # exe_wfc_real: complex array (nstates,nspinor_electron,nspinor_hole, FFTx, FFTy, FFTz). - # Please note that you must take absoulte square and contract the spinor dimension to the density. + Returns: + tuple: (supercell_latvecs, atom_nums, atom_pos, exe_wfc_real) + Natoms = (natom_in_unit_call * Nsupercell) + 1 (+1 due to hole/electron) + - supercell_latvecs (numpy.ndarray): Supercell lattice vectors [3,3] + - atom_nums (numpy.ndarray): Atomic numbers [Natoms] + - atom_pos (numpy.ndarray): Atomic positions in cartisian units [Natoms,3] + - exe_wfc_real (numpy.ndarray): Wavefunction in real space + [nstates, nspinor_electron, nspinor_hole, FFTx, FFTy, FFTz] + """ + # # if block_size < 1: print('Warning: Wrong block_size. setting to 1') @@ -74,7 +183,16 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, hole_bnds = [bse_bnds[0],bse_bnds[0]+nv] elec_bnds = [bse_bnds[0]+nv,bse_bnds[1]] - + # + if ares: + # if anti-resonant part, we swap the electron and hole bands + # first transpose c,v dimensions + Akcv = Akcv.transpose(0,1,3,2) + tmp = hole_bnds + hole_bnds = elec_bnds + elec_bnds = tmp + nstates, nk, nc, nv = Akcv.shape + lat_vec = wfcdb.ydb.lat.T blat = np.linalg.inv(lat_vec) gvecs_iBZ_idx = [] @@ -98,13 +216,23 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, fft_box[i] = max([fft_box[i], max_fft_idx[i] - min_fft_idx[i] + 3]) # Compute nstates, nk, Nx, Ny, Nz object - print("FFT Box : ",fft_box[0], fft_box[1], fft_box[2]) + if out_res is None : print("Wfc FFT Grid : ",fft_box[0], fft_box[1], fft_box[2]) # ktree = build_ktree(wfcdb.kBZ) # nspinorr = wfcdb.nspinor - exe_wfc_real = np.zeros((nstates, nspinorr, nspinorr, np.prod(supercell), - fft_box[0], fft_box[1], fft_box[2]),dtype=np.complex64) + if ares and out_res is not None: + exe_wfc_real = out_res.reshape(nstates, nspinorr, nspinorr, + supercell[0],fft_box[0], + supercell[1],fft_box[1], + supercell[2],fft_box[2]) + else: + exe_wfc_real = np.zeros((nstates, nspinorr, nspinorr, + supercell[0],fft_box[0], + supercell[1],fft_box[1], + supercell[2],fft_box[2]), + dtype=np.complex64) + # Lx = np.arange(supercell[0],dtype=int) Ly = np.arange(supercell[1],dtype=int) Lz = np.arange(supercell[2],dtype=int) @@ -130,6 +258,8 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, nrem = nk%block_size if nrem > 0: nblks = nblks+1 # + pbar = tqdm(total=nk, desc="Ex-wf") + # for ibk in range(nblks): ikstart = ibk*block_size ikstop = min(ikstart + block_size,nk) @@ -209,6 +339,10 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] exp_tmp_kL[ik-ikstart] = np.exp(1j*2*np.pi*np.einsum('...x,x->...',Lsupercells,ft_kvec)) + # + # update progess bar + pbar.update(1) + # ## perform gemm operation total_gemms_t = nstates*nspinorr**2 exp_tmp_kL_tmp = exp_tmp_kL.reshape(len(exp_tmp_kL),-1)[:(ikstop-ikstart)].T @@ -217,19 +351,17 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # for igemms in range(total_gemms_t): ii, jj, kk = np.unravel_index(igemms, (nstates,nspinorr,nspinorr)) - # NM : It is not nice to create an large temporary array again. but numpy does support - # C += A@B call like that blas has. - exe_wfc_real[ii, jj, kk ] += (exp_tmp_kL_tmp @ exe_tmp_wf_tmp[ii, jj, kk ]).reshape( - np.prod(supercell),fft_box[0], fft_box[1], fft_box[2]) - - exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr**2,supercell[0],supercell[1],supercell[2], - fft_box[0], fft_box[1], fft_box[2]) + Ctmp = (exp_tmp_kL_tmp @ exe_tmp_wf_tmp[ii, jj, kk ]) + Ctmp = Ctmp.reshape(supercell[0], supercell[1], supercell[2], + fft_box[0], fft_box[1], fft_box[2]) + exe_wfc_real[ii, jj, kk ] += Ctmp.transpose(0,3,1,4,2,5) # - exe_wfc_real = exe_wfc_real.transpose(0,1,2,5,3,6,4,7).reshape(nstates,nspinorr,nspinorr, - supercell[0]*fft_box[0], supercell[1]*fft_box[1], - supercell[2]*fft_box[2]) + exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr,nspinorr, + supercell[0]*fft_box[0], + supercell[1]*fft_box[1], + supercell[2]*fft_box[2]) exe_wfc_real *= (1.0/np.prod(supercell)) - + # # compute postioon of fixed particle in cart units fixed_postion_cc = lat_vec@fixed_postion Lsupercells = Lsupercells.reshape(-1,3)#/np.array(supercell)[None,:] @@ -245,56 +377,56 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # return supercell_latvecs,atom_nums,atom_pos,exe_wfc_real -def compute_exc_wfc_real(path='.', bse_dir='SAVE', iqpt=1, nstates=[1], - fixed_postion=[0,0,0], fix_particle='h', aveg=True, supercell=[1,1,1], - wfcCutoffRy=-1, phase=False, block_size=256): - # - lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) - filename = 'ndb.BS_diago_Q%d' % (iqpt) - excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, - folder=os.path.join(path, bse_dir), - Load_WF=True, neigs=max(nstates)) - # Load the wavefunction database - wfdb = YamboWFDB(path=path, latdb=lattice, - bands_range=[np.min(excdb.table[:, 1]) - 1, - np.max(excdb.table[:, 2])]) - # - Akcv = excdb.get_Akcv()[min(nstates)-1:max(nstates)] - excQpt = excdb.car_qpoint - # - # Convert the q-point to crystal coordinates - Qpt = wfdb.ydb.lat @ excQpt +#def compute_exc_wfc_real(path='.', bse_dir='SAVE', iqpt=1, nstates=[1], +# fixed_postion=[0,0,0], fix_particle='h', aveg=True, supercell=[1,1,1], +# wfcCutoffRy=-1, phase=False, block_size=256): +# # +# lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1')) +# filename = 'ndb.BS_diago_Q%d' % (iqpt) +# excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, +# folder=os.path.join(path, bse_dir), +# Load_WF=True, neigs=max(nstates)) +# # Load the wavefunction database +# wfdb = YamboWFDB(path=path, latdb=lattice, +# bands_range=[np.min(excdb.table[:, 1]) - 1, +# np.max(excdb.table[:, 2])]) +# # +# Akcv = excdb.get_Akcv()[min(nstates)-1:max(nstates)] +# excQpt = excdb.car_qpoint +# # +# # Convert the q-point to crystal coordinates +# Qpt = wfdb.ydb.lat @ excQpt - sc_latvecs, atom_nums, atom_pos, real_wfc = ex_wf2Real(Akcv, Qpt, wfdb, [np.min(excdb.table[:, 1]), - np.max(excdb.table[:, 2])], fixed_postion=fixed_postion, - fix_particle=fix_particle, supercell=supercell, - wfcCutoffRy=wfcCutoffRy, block_size=block_size) - # - # - nstates_range = np.arange(nstates[0],nstates[1],dtype=int) - density = np.abs(real_wfc)**2 - - if fix_particle == 'h': name_file = 'electron' - else: name_file = 'hole' - if real_wfc.shape[1] != 1: - print("phase plot only works for nspin = 1 and nspinor == 1") - phase = False - if phase: - phase = np.sign(real_wfc.real) #np.sign(np.angle(real_wfc)) - density *= phase - if aveg: - real_wfc = np.sum(density,axis=(0,1,2)) - real_wfc = real_wfc/np.max(np.abs(real_wfc)) - write_cube('exe_wf_avg_%s_%d-%d.cube' %(name_file,nstates[0],nstates[1]), - real_wfc, sc_latvecs, atom_pos, atom_nums, header='Real space exciton wavefunction') - else: - real_wfc = np.sum(density,axis=(1,2)) - for i in range(len(real_wfc)): - real_wfc1 = real_wfc[i]/np.max(np.abs(real_wfc[i])) - write_cube('exe_wf_%s_%d.cube' %(name_file,nstates[i]), real_wfc1, sc_latvecs, - atom_pos, atom_nums, header='Real space exciton wavefunction') +# sc_latvecs, atom_nums, atom_pos, real_wfc = ex_wf2Real(Akcv, Qpt, wfdb, [np.min(excdb.table[:, 1]), +# np.max(excdb.table[:, 2])], fixed_postion=fixed_postion, +# fix_particle=fix_particle, supercell=supercell, +# wfcCutoffRy=wfcCutoffRy, block_size=block_size) +# # +# # +# nstates_range = np.arange(nstates[0],nstates[1],dtype=int) +# density = np.abs(real_wfc)**2 + +# if fix_particle == 'h': name_file = 'electron' +# else: name_file = 'hole' +# if real_wfc.shape[1] != 1: +# print("phase plot only works for nspin = 1 and nspinor == 1") +# phase = False +# if phase: +# phase = np.sign(real_wfc.real) #np.sign(np.angle(real_wfc)) +# density *= phase +# if aveg: +# real_wfc = np.sum(density,axis=(0,1,2)) +# real_wfc = real_wfc/np.max(np.abs(real_wfc)) +# write_cube('exe_wf_avg_%s_%d-%d.cube' %(name_file,nstates[0],nstates[1]), +# real_wfc, sc_latvecs, atom_pos, atom_nums, header='Real space exciton wavefunction') +# else: +# real_wfc = np.sum(density,axis=(1,2)) +# for i in range(len(real_wfc)): +# real_wfc1 = real_wfc[i]/np.max(np.abs(real_wfc[i])) +# write_cube('exe_wf_%s_%d.cube' %(name_file,nstates[i]), real_wfc1, sc_latvecs, +# atom_pos, atom_nums, header='Real space exciton wavefunction') -#if __name__ == "__main__": +##if __name__ == "__main__": diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index a5e12482..9b6eb27f 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -24,6 +24,8 @@ from yambopy.dbs.latticedb import YamboLatticeDB from yambopy.dbs.electronsdb import YamboElectronsDB from yambopy.dbs.qpdb import YamboQPDB +from yambopy.io.cubetools import write_cube +from yambopy.bse.realSpace_excitonwf import ex_wf2Real class ExcitonList(): """ @@ -278,6 +280,70 @@ def get_Akcv(self): # self.Akcv = eig_wfcs_returned return self.Akcv + + def real_wf_to_cube(iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], degen_tol=1e-2, + wfcCutoffRy=-1, fix_particle='h', phase=False, block_size=256): + """ + Function to compute and save real-space exciton wavefunctions and + dump to cube file + + Args: + iexe: index of excitonic states (python indexing. so starts with 0) + wfcb: wavefunction database. + fixed_postion (list): Position of fixed particle in crystal coordinates + supercell (list): Supercell dimensions [nx,ny,nz] + degen_tol (float): degeneracy threshold (in eV). default 0.01 eV + fix_particle (str): 'e' to fix electron, 'h' to fix hole (default) + wfcCutoffRy (float): Wavefunction cutoff in Rydberg (-1 for full cutoff) + phase (bool): If True, include phase information i.e multiply the density with + sign of real part of the wavefunction + block_size (int): Block size for memory-efficient computation. leave it to default + unless you are in exteremely low memory situation. + + Returns: + None (write cube file to disk) + """ + # + # first get all degenerate states + iexe_degen_states = np.array(self.get_degenerate(iexe+1,eps=degen_tol))-1 + # nicely arrange eigvectors to Akcv + Akcv = self.get_Akcv()[iexe_degen_states] + excQpt = self.car_qpoint + # Convert the q-point to crystal coordinates + Qpt = wfdb.ydb.lat @ excQpt + # + if fix_particle == 'h': name_file = 'electron' + else: name_file = 'hole' + + if phase and real_wfc.shape[1] != 1: + print("phase plot only works for nspin = 1 and nspinor == 1") + phase = False + if phase and len(iexe_degen_states) > 1: + phase = False + print("Warning: phase plots donot work for degenerate states") + + print('Computing exciton wavefunction (%s density) to real space.' %(name_file)) + sc_latvecs, atom_nums, atom_pos, real_wfc = ex_wf2Real(Akcv, Qpt, wfdb, [np.min(self.table[:, 1]), + np.max(self.table[:, 2])], fixed_postion=fixed_postion, + fix_particle=fix_particle, supercell=supercell, + wfcCutoffRy=wfcCutoffRy, block_size=block_size) + # Compute the absoulte^2 + density = np.abs(real_wfc)**2 + # Multiply with phase if necessary + if phase: + phase = np.sign(real_wfc.real) #np.sign(np.angle(real_wfc)) + density *= phase + # + # sum over spinor indices and degenerate states + real_wfc = np.sum(density,axis=(0,1,2)) + # normalize with max value + real_wfc = real_wfc/np.max(np.abs(real_wfc)) + # write to cube file + print('Writing to .cube file') + write_cube('exe_wf_%s_%d.cube' %(name_file,iexe+1), + real_wfc, sc_latvecs, atom_pos, atom_nums, + header='Real space exciton wavefunction') + def get_nondegenerate(self,eps=1e-4): """ @@ -325,12 +391,9 @@ def get_degenerate(self,index,eps=1e-4): Args: eps: maximum energy difference to consider the two excitons degenerate in eV """ - energy = self.eigenvalues[index-1] - excitons = [] - for n,e in enumerate(self.eigenvalues): - if np.isclose(energy,e,atol=eps): - excitons.append(n+1) - return excitons + energy = self.eigenvalues[index-1].real + excitons = np.where(np.isclose(self.eigenvalues.real, energy, atol=eps))[0] + 1 + return excitons.tolist() def exciton_bs(self,energies,path,excitons=(0,),debug=False): """ From de45ca2eaa268d999d3b29412d638afff46e45ef Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 13 Apr 2025 23:42:45 +0200 Subject: [PATCH 44/73] tiny fixes --- yambopy/bse/realSpace_excitonwf.py | 10 ++++------ yambopy/dbs/excitondb.py | 3 ++- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 31d9d202..f69c2989 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -1,9 +1,6 @@ ### Compute real space exction wavefunction when hole/electron is fixed. import numpy as np from yambopy.kpoints import build_ktree, find_kpt -from yambopy.dbs.excitondb import YamboExcitonDB -from yambopy.dbs.latticedb import YamboLatticeDB -from yambopy.dbs.wfdb import YamboWFDB from tqdm import tqdm import os @@ -13,6 +10,7 @@ from yambopy.dbs.excitondb import YamboExcitonDB from yambopy.dbs.latticedb import YamboLatticeDB from yambopy.dbs.wfdb import YamboWFDB +import os iqpt = 1 # qpt index of exciton @@ -25,11 +23,11 @@ # filename = 'ndb.BS_diago_Q%d' % (iqpt) excdb = YamboExcitonDB.from_db_file(lattice, filename=filename, - folder=os.path.join(path, bse_dir) + folder=os.path.join('.', 'GW_BSE'), neigs = 20) #Load the wavefunction database -wfdb = YamboWFDB(path=path, latdb=lattice, +wfdb = YamboWFDB(path='.', latdb=lattice, bands_range=[np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])]) @@ -38,7 +36,7 @@ # 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, fixed_postion=[0,0,0], supercell=[1,1,1], +excdb.real_wf_to_cube(iexe=2, wfdb=wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], degen_tol=0.01, wfcCutoffRy=80, fix_particle='h') ## .cube will be dumped and use vesta to visualize it ! diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index 9b6eb27f..3ab647f0 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -281,7 +281,7 @@ def get_Akcv(self): self.Akcv = eig_wfcs_returned return self.Akcv - def real_wf_to_cube(iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], degen_tol=1e-2, + def real_wf_to_cube(self, iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], degen_tol=1e-2, wfcCutoffRy=-1, fix_particle='h', phase=False, block_size=256): """ Function to compute and save real-space exciton wavefunctions and @@ -306,6 +306,7 @@ def real_wf_to_cube(iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], degen_ # # first get all degenerate states iexe_degen_states = np.array(self.get_degenerate(iexe+1,eps=degen_tol))-1 + print("Degenerate states: ",iexe_degen_states+1) # nicely arrange eigvectors to Akcv Akcv = self.get_Akcv()[iexe_degen_states] excQpt = self.car_qpoint From 15b59a91874271884d19bda14f64eb49b2bb0094 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 15 Apr 2025 08:27:23 +0200 Subject: [PATCH 45/73] add kpatch function --- yambopy/bse/realSpace_excitonwf.py | 4 ++-- yambopy/dbs/excitondb.py | 4 +++- yambopy/kpoints.py | 35 ++++++++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 3 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index f69c2989..8dbc8dcb 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -31,7 +31,7 @@ bands_range=[np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])]) -## plot the exciton wavefunction with hole fixed at [0,0,0] +## plot the exciton wavefunction with hole fixed at [0,0,0] (crystal coordinates) # in a [1,1,1] supercell with 80 Ry wf 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 ) @@ -39,6 +39,7 @@ excdb.real_wf_to_cube(iexe=2, wfdb=wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], degen_tol=0.01, wfcCutoffRy=80, fix_particle='h') +## fix_particle = 'e' if you want to fix electron and plot hole density ## .cube will be dumped and use vesta to visualize it ! """ @@ -316,7 +317,6 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, fx_wfc = wfc_elec fx_gvec = gvecs_elec # compute - ## Fix compute the fixed particle wfc in real space. ## NM : Donot perform FFT as we only need it for one point. exp_fx = np.exp(2*np.pi*1j*((fx_gvec + fx_kvec[None,:])@fixed_postion)) fx_wfc *= exp_fx[None,None,None,:] diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index 3ab647f0..be22d8d0 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -338,7 +338,9 @@ def real_wf_to_cube(self, iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], # sum over spinor indices and degenerate states real_wfc = np.sum(density,axis=(0,1,2)) # normalize with max value - real_wfc = real_wfc/np.max(np.abs(real_wfc)) + max_normalize_val = np.max(np.abs(real_wfc)) + print('Max Normalization value: ',max_normalize_val) + real_wfc *= (1.0/max_normalize_val) # write to cube file print('Writing to .cube file') write_cube('exe_wf_%s_%d.cube' %(name_file,iexe+1), diff --git a/yambopy/kpoints.py b/yambopy/kpoints.py index 2247968c..e01cdf1e 100644 --- a/yambopy/kpoints.py +++ b/yambopy/kpoints.py @@ -213,4 +213,39 @@ def find_kpt(tree, kpt_search, tol=1e-5): return idx # Return the index of the found k-point +def find_kpatch(kpts, kcentre, kdist, lat_vecs): + """ + find set of kpoints around the kcentre with in kdist + + Parameters + ---------- + kpts : kpoints in crystal coordinates (nk,3) + kcentre : kpoint centre in crystal coordinates (3) + kdist : distance around kcentre to be considered in atomic units + i.e 1/bohr. + lat_vecs: lattice vectors. ith lattice vector is ai = a[:,i] + Returns + ------- + int array + Indices of kpoints in kpts array which satify the given condition i.e + | k - kcentre + G0| <= kdist, where G0 is reciprocal lattice vector to bring to BZ + """ + # + blat = 2*np.pi*np.linalg.inv(lat_vecs) + kdiff = kpts-kcentre[None,:] + kdiff = kdiff-np.floor(kdiff) + # + tmp_arr = np.array([-3, -2, -1, 0, 1, 2, 3]) + nG0 = len(tmp_arr) + G0 = np.zeros((nG0,nG0,nG0,3)) + G0[...,0], G0[...,1], G0[...,2] = np.meshgrid(tmp_arr, tmp_arr, + tmp_arr, indexing='ij') + G0 = G0.reshape(-1,3) + kdiff = kdiff[:,None,:]-G0[None,:,:] + kdiff = kdiff.reshape(-1,3)@blat + kdiff = np.linalg.norm(kdiff,axis=-1).reshape(len(kpts),-1) + kdiff = np.min(kdiff,axis=-1) + return np.where(kdiff <= kdist)[0] + + From 85af27ad9746f5859e0bf6d16db23cefb46045e6 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 17 Apr 2025 15:18:04 +0200 Subject: [PATCH 46/73] . --- yambopy/bse/realSpace_excitonwf.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 8dbc8dcb..6256c508 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -129,8 +129,9 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, the default is 256 which is very good but uses more memory. # decrease it when you run into memory issues ares (bool): If True, compute anti-resonant part - out_res (numpy.ndarray): If provided and ares=True, adds to this array + out_res (numpy.ndarray): Adds to this array and is returned instead of internally creating. + Make sure it its consistant. (no internal checking done) Returns: tuple: (supercell_latvecs, atom_nums, atom_pos, exe_wfc_real) @@ -220,7 +221,7 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, ktree = build_ktree(wfcdb.kBZ) # nspinorr = wfcdb.nspinor - if ares and out_res is not None: + if out_res is not None: exe_wfc_real = out_res.reshape(nstates, nspinorr, nspinorr, supercell[0],fft_box[0], supercell[1],fft_box[1], From 91b97f4fcbdbfca30157513d4c0f6d61159df34c Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 17 Apr 2025 20:02:00 +0200 Subject: [PATCH 47/73] remove [0,1) creteria --- yambopy/bse/realSpace_excitonwf.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 6256c508..3d7c75fa 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -173,10 +173,11 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, assert nk == nkBZ, "kpoint mismatch" # - ## Bring the positon of hole/electron in the centre unit cell, i.e make in between [0,1) - fixed_postion = np.array(fixed_postion) - np.floor(fixed_postion) - fixed_postion = (fixed_postion + 1e-6)%1 - fixed_postion = fixed_postion - 1e-6 + # ## Bring the positon of hole/electron in the centre unit cell, i.e make in between [0,1) + # fixed_postion = np.array(fixed_postion) - np.floor(fixed_postion) + # fixed_postion = (fixed_postion + 1e-6)%1 + # fixed_postion = fixed_postion - 1e-6 + fixed_postion = np.array(fixed_postion) fixed_postion += np.array(supercell)//2 From cbf54644b393ea7ecfe6641fabf21d1ed739eaef Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 22 May 2025 20:58:47 +0200 Subject: [PATCH 48/73] update progess --- yambopy/bse/exciton_matrix_elements.py | 2 ++ yambopy/bse/realSpace_excitonwf.py | 3 ++ yambopy/bse/rotate_excitonwf.py | 2 ++ yambopy/dbs/wfdb.py | 11 +++++-- yambopy/tools/function_profiler.py | 40 ++++++++++++++++++++++++++ 5 files changed, 55 insertions(+), 3 deletions(-) create mode 100644 yambopy/tools/function_profiler.py diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index 38cb7ffd..eed233c8 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -8,7 +8,9 @@ ### 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. diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 3d7c75fa..db7ef3ee 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -1,6 +1,7 @@ ### Compute real space exction wavefunction when hole/electron is fixed. import numpy as np from yambopy.kpoints import build_ktree, find_kpt +from yambopy.tools.function_profiler import func_profile from tqdm import tqdm import os @@ -45,6 +46,7 @@ ## +@func_profile def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, fix_particle='h', supercell=[1,1,1], wfcCutoffRy=-1, block_size=256): @@ -102,6 +104,7 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, +@func_profile def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, fix_particle='h', supercell=[1,1,1], wfcCutoffRy=-1, block_size=256, ares=False, out_res=None): diff --git a/yambopy/bse/rotate_excitonwf.py b/yambopy/bse/rotate_excitonwf.py index 72a2ad8e..7595aa97 100644 --- a/yambopy/bse/rotate_excitonwf.py +++ b/yambopy/bse/rotate_excitonwf.py @@ -1,6 +1,8 @@ import numpy as np from yambopy.kpoints import build_ktree, find_kpt +from yambopy.tools.function_profiler import func_profile +@func_profile def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=None): """ Rotate the exciton wavefunction Ak using symmetry operations. diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 4b55c723..c354db77 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -25,6 +25,7 @@ except ImportError as e: from scipy.spatial import KDTree from yambopy.kpoints import build_ktree, find_kpt +from yambopy.tools.function_profiler import func_profile class YamboWFDB: """ @@ -389,6 +390,7 @@ def rotate_wfc(self, ik, isym): time_rev = (isym >= len(self.ydb.sym_car) / (1 + int(np.rint(self.ydb.time_rev)))) return self.apply_symm(kvec, wfc_k, gvecs_k, time_rev, sym_mat) + @func_profile def apply_symm(self, kvec, wfc_k, gvecs_k, time_rev, sym_mat, frac_vec=np.array([0, 0, 0])): """ Apply symmetry to wavefunctions. @@ -426,6 +428,7 @@ def apply_symm(self, kvec, wfc_k, gvecs_k, time_rev, sym_mat, frac_vec=np.array( return [Rkvec, wfc_rot, gvec_rot] + @func_profile def to_real_space(self, wfc_tmp, gvec_tmp, grid=[]): """ Convert wavefunctions from G-space to real space. @@ -516,6 +519,7 @@ def get_BZ_wf(self, ik): # return [self.wf_bz[ik][..., :self.ngBZ[ik]], self.g_bz[ik, :self.ngBZ[ik], :]] + @func_profile def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): """ Computes the symmetry-adapted matrix elements < Rk | U(R) | k >. @@ -601,7 +605,8 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): # return Dmat -def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gtree=-1): +@func_profile +def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gtree=None): """ Computes the inner product between two wavefunctions in reciprocal space. @@ -620,7 +625,7 @@ def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gt gvec_ket : ndarray Miller indices of the ket wavefunction (ng, 3) in reduced coordinates. ket_Gtree : scipy.spatial._kdtree.KDTree (optional) - Kdtree for gvec_ket. leave it or give -1 to internally build one + Kdtree for gvec_ket. leave it or give None to internally build one # Returns ------- @@ -640,7 +645,7 @@ def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gt if np.max(np.abs(kdiff)) > 1e-5: return np.zeros((nspin, nbnd, nbnd),dtype=wfc_ket.dtype) # Construct KDTree for nearest-neighbor search in G-vectors - if type(ket_Gtree) != scipy.spatial._kdtree.KDTree: + if ket_Gtree is None: ket_Gtree = KDTree(gvec_ket) gbra_shift = gvec_bra + G0[None,:] ## get the nearest indices and their distance diff --git a/yambopy/tools/function_profiler.py b/yambopy/tools/function_profiler.py new file mode 100644 index 00000000..a43d2038 --- /dev/null +++ b/yambopy/tools/function_profiler.py @@ -0,0 +1,40 @@ +import time +from functools import wraps + +# decarator to measure number of calls and wall time in seconds of function +""" +@func_profile +def abc(xxx): + np.log(np.sin(xxx)+1j+np.cos(xxx)+xxx**2) + +x = np.linspace(0,10,1000000) + +for i in range(10): + abc(x) + +# number of calls to abc +print(abc.call_count) + +# total time spent in abc +print(abc.total_time) + +""" +def func_profile(func): + @wraps(func) + def wrapper(*args, **kwargs): + start_time = time.perf_counter() + + # Increment call count (initialize if first call) + if not hasattr(wrapper, "call_count"): + wrapper.call_count = 0 + wrapper.call_count += 1 + + result = func(*args, **kwargs) # Call original function + + # Accumulate time (initialize if first call) + if not hasattr(wrapper, "total_time"): + wrapper.total_time = 0.0 + wrapper.total_time += time.perf_counter() - start_time + + return result + return wrapper From 5b49f461e828c0793f75c231d33bc95613a44665 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 29 May 2025 14:29:30 +0200 Subject: [PATCH 49/73] update wfdb copyright --- yambopy/dbs/wfdb.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index e6f40e90..3f04bced 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -1,8 +1,6 @@ -# Copyright (c) 2018, Henrique Miranda +# Copyright (c) 2025, Muralidhar Nalabothula # All rights reserved. # -# This file is part of the yambopy project -# # Author: MN from yambopy import * From 50d34bfb92dcdd8fb41fa79033a94e0222b0c4ae Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 2 Jun 2025 18:15:46 +0200 Subject: [PATCH 50/73] more roboust nc_char to string --- yambopy/letzelphc_interface/lelphcdb.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/yambopy/letzelphc_interface/lelphcdb.py b/yambopy/letzelphc_interface/lelphcdb.py index f8463e90..a0150b15 100644 --- a/yambopy/letzelphc_interface/lelphcdb.py +++ b/yambopy/letzelphc_interface/lelphcdb.py @@ -50,9 +50,15 @@ def __init__(self,filename,read_all=True,div_by_energies=True,verbose=False): self.nsym = database.dimensions['nsym_ph'].size conv = database['convention'][...].data - convlist = [iconv.decode('utf-8') for iconv in conv] - conv = '' - for iconv in convlist: conv = conv + iconv + if isinstance(conv, np.ndarray): + if conv.dtype.kind == 'S': # Byte strings (C chars) + conv = conv.tobytes().decode('utf-8').strip() + else: + conv = str(conv) # Fallback for non-string arrays + elif isinstance(conv, bytes): + conv = conv.decode('utf-8').strip() + else: + conv = str(conv).strip() stard_conv = False if conv.strip() == 'standard': print("Convention used in Letzelphc : k -> k+q (standard)") From 086bade308fc1ec3e54051046d42231b18241ed3 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 8 Jul 2025 13:21:18 +0200 Subject: [PATCH 51/73] - --- yambopy/letzelphc_interface/lelphcdb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/letzelphc_interface/lelphcdb.py b/yambopy/letzelphc_interface/lelphcdb.py index a0150b15..92f9d34f 100644 --- a/yambopy/letzelphc_interface/lelphcdb.py +++ b/yambopy/letzelphc_interface/lelphcdb.py @@ -8,7 +8,7 @@ class LetzElphElectronPhononDB(): """ Python class to read the electron-phonon matrix elements from LetzElPhC. - About LetzElPhC: https://github.com/yambo-code/LetzElPhC/tree/main + About LetzElPhC: https://gitlab.com/lumen-code/LetzElPhC By default it reads the full database g(q,k,m,s,b1,b2) including phonon energies. From c7e47f918f5c14cac99810c7bd948739f35c41cc Mon Sep 17 00:00:00 2001 From: Murali Date: Sat, 12 Jul 2025 12:32:15 +0200 Subject: [PATCH 52/73] a --- yambopy/io/cubetools.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/yambopy/io/cubetools.py b/yambopy/io/cubetools.py index 934619e2..d5f11100 100644 --- a/yambopy/io/cubetools.py +++ b/yambopy/io/cubetools.py @@ -11,10 +11,10 @@ def write_cube(filename, data, lat_vec, atom_pos, atomic_num, origin=np.zeros(3) atomic_num: atomic numbers of elements """ ## atomic positions must be in [0,1) - apos_crys = atom_pos@np.linalg.inv(lat_vec.T) + apos_crys = atom_pos@lat_vec.T apos_crys = apos_crys-np.floor(apos_crys) apos_crys = (apos_crys+1e-6)%1 - apos_new = apos_crys@lat_vec.T + apos_new = apos_crys@np.linalg.inv(lat_vec.T) with open(filename, "w") as cube: cube.write("# Cube file generated using YamboPy\n") @@ -35,5 +35,4 @@ def write_cube(filename, data, lat_vec, atom_pos, atomic_num, origin=np.zeros(3) for iz in range(data.shape[2]): cube.write("%.6f " %(data[ix,iy,iz])) if iz % 6 == 5 : cube.write("\n") - cube.write("\n") From 9a6351c52ed9e8f19fee8b9bd43f8a106eac66c9 Mon Sep 17 00:00:00 2001 From: Murali Date: Sat, 12 Jul 2025 12:33:49 +0200 Subject: [PATCH 53/73] revert --- yambopy/io/cubetools.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/yambopy/io/cubetools.py b/yambopy/io/cubetools.py index d5f11100..934619e2 100644 --- a/yambopy/io/cubetools.py +++ b/yambopy/io/cubetools.py @@ -11,10 +11,10 @@ def write_cube(filename, data, lat_vec, atom_pos, atomic_num, origin=np.zeros(3) atomic_num: atomic numbers of elements """ ## atomic positions must be in [0,1) - apos_crys = atom_pos@lat_vec.T + apos_crys = atom_pos@np.linalg.inv(lat_vec.T) apos_crys = apos_crys-np.floor(apos_crys) apos_crys = (apos_crys+1e-6)%1 - apos_new = apos_crys@np.linalg.inv(lat_vec.T) + apos_new = apos_crys@lat_vec.T with open(filename, "w") as cube: cube.write("# Cube file generated using YamboPy\n") @@ -35,4 +35,5 @@ def write_cube(filename, data, lat_vec, atom_pos, atomic_num, origin=np.zeros(3) for iz in range(data.shape[2]): cube.write("%.6f " %(data[ix,iy,iz])) if iz % 6 == 5 : cube.write("\n") + cube.write("\n") From cb3f9b5d688857a22acd7349e3d6d35b692c6f26 Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 13 Jul 2025 19:32:46 +0200 Subject: [PATCH 54/73] drr wip --- yambopy/bse/realSpace_excitonwf.py | 31 ++++++++++++++++++------------ yambopy/dbs/wfdb.py | 13 ++++++++++++- 2 files changed, 31 insertions(+), 13 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index db7ef3ee..7731d44d 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -175,16 +175,8 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, assert nc + nv == bse_bnds[1]-bse_bnds[0], "Band mismatch" assert nk == nkBZ, "kpoint mismatch" # - - # ## Bring the positon of hole/electron in the centre unit cell, i.e make in between [0,1) - # fixed_postion = np.array(fixed_postion) - np.floor(fixed_postion) - # fixed_postion = (fixed_postion + 1e-6)%1 - # fixed_postion = fixed_postion - 1e-6 fixed_postion = np.array(fixed_postion) - - fixed_postion += np.array(supercell)//2 - - + # hole_bnds = [bse_bnds[0],bse_bnds[0]+nv] elec_bnds = [bse_bnds[0]+nv,bse_bnds[1]] # @@ -222,7 +214,21 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # Compute nstates, nk, Nx, Ny, Nz object if out_res is None : print("Wfc FFT Grid : ",fft_box[0], fft_box[1], fft_box[2]) # - ktree = build_ktree(wfcdb.kBZ) + ## + # find the nearest fft grid point. + fx_pnt_int = np.floor(fixed_postion) + fixed_postion -= fx_pnt_int + fixed_postion = np.round(fixed_postion * fft_box) / fft_box + fixed_postion += fx_pnt_int + # shift the position of hole to middle of supercell + fixed_postion += np.array(supercell)//2 + # + if fix_particle == 'h': + print("Position of the hole is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) + if fix_particle == 'e': + print("Position of the electron is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) + # + ktree = wfcdb.ktree #build_ktree(wfcdb.kBZ) # nspinorr = wfcdb.nspinor if out_res is not None: @@ -331,7 +337,8 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, #if ft_ikpt not in prev_ikpts: ft_wfcr = wfcdb.to_real_space(ft_wfc.reshape(-1,nspinorr,ng),ft_gvec, grid=fft_box) ft_wfcr = ft_wfcr.reshape(ns1,nbndc,nspinorr,fft_box[0],fft_box[1],fft_box[2]) - exp_kx_r = np.exp(2*np.pi*1j*FFFboxs.reshape(-1,3)@ft_kvec) + exp_kx_r = np.exp(2*np.pi*1j*FFFboxs.reshape(-1,3)@ft_kvec).reshape(FFFboxs.shape[:3]) + ft_wfcr *= exp_kx_r[None,None,None,...] # if fix_particle == 'h': np.einsum('ncv,vy,cxijk->nxyijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], @@ -340,7 +347,7 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, np.einsum('ncv,cx,vyijk->nxyijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], optimize=True,out=exe_tmp_wf[:,:,:,ik-ikstart]) # - exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] + #exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] exp_tmp_kL[ik-ikstart] = np.exp(1j*2*np.pi*np.einsum('...x,x->...',Lsupercells,ft_kvec)) # # update progess bar diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 3f04bced..50fdea03 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -128,6 +128,7 @@ def read(self, bands_range=[], latdb=None): # K-points in BZ (crystal units) self.kBZ = self.ydb.iku_kpoints / lat_param[None, :] self.kBZ = self.kBZ @ lat_vec + self.ktree = build_ktree(self.kBZ) # G-vectors in cartesian units G_vec = ns_db1['G-VECTORS'][...].data.T @@ -214,6 +215,13 @@ def assert_k_inrange(self, ik): def assert_bnd_range(self, ib): """Assert that the band index is valid.""" assert 0 <= ib < self.nbands, "Invalid band index" + + def kptBZidx(self, kpts): + """ + return the index of the kpoint or kpoints in the wavefunction db. + The kpts must be in crystal (reduced) coordinates. + """ + return find_kpt(self.ktree, kpts) def get_spin_projections(self, ik, ib=-1, s_z=np.array([[1, 0], [0, -1]])): """ @@ -493,6 +501,9 @@ def expand_fullBZ(self): self.kBZ[i] = kbz self.wf_bz[i][...,:ng_t] = w_t self.g_bz[i][:ng_t,:] = g_t + # + self.ktree = build_ktree(self.kBZ) + return def get_BZ_kpt(self, ik): """ @@ -563,7 +574,7 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): frac_vec = np.zeros((len(symm_mat),3),dtype=symm_mat.dtype) time_rev = int(np.rint(self.ydb.time_rev)) # - ktree = build_ktree(self.kBZ) + ktree = self.ktree #build_ktree(self.kBZ) Dmat = [] nsym = len(symm_mat) kpt_idx = self.ydb.kpoints_indexes From c73eb9038e20b471af7dce085785b60d3917518a Mon Sep 17 00:00:00 2001 From: Murali Date: Sun, 13 Jul 2025 20:02:04 +0200 Subject: [PATCH 55/73] drr wip --- yambopy/bse/realSpace_excitonwf.py | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 7731d44d..8bb59aa0 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -176,6 +176,15 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, assert nk == nkBZ, "kpoint mismatch" # fixed_postion = np.array(fixed_postion) + fixed_postion += np.array(supercell)//2 + # + if fix_particle == 'h': + print("Position of the hole is set to (reduced units) : ", + fixed_postion[0], fixed_postion[1], fixed_postion[2]) + elif fix_particle == 'e': + print("Position of the electron is set to (reduced units) : ", + fixed_postion[0], fixed_postion[1], fixed_postion[2]) + # # hole_bnds = [bse_bnds[0],bse_bnds[0]+nv] elec_bnds = [bse_bnds[0]+nv,bse_bnds[1]] @@ -214,21 +223,7 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # Compute nstates, nk, Nx, Ny, Nz object if out_res is None : print("Wfc FFT Grid : ",fft_box[0], fft_box[1], fft_box[2]) # - ## - # find the nearest fft grid point. - fx_pnt_int = np.floor(fixed_postion) - fixed_postion -= fx_pnt_int - fixed_postion = np.round(fixed_postion * fft_box) / fft_box - fixed_postion += fx_pnt_int - # shift the position of hole to middle of supercell - fixed_postion += np.array(supercell)//2 - # - if fix_particle == 'h': - print("Position of the hole is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) - if fix_particle == 'e': - print("Position of the electron is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) - # - ktree = wfcdb.ktree #build_ktree(wfcdb.kBZ) + ktree = wfcdb.ktree # nspinorr = wfcdb.nspinor if out_res is not None: From 0f55ca6ad1287451bc0a5e28306588e18ec4ba4e Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 14 Jul 2025 12:21:59 +0200 Subject: [PATCH 56/73] Revert "drr wip" This reverts commit c73eb9038e20b471af7dce085785b60d3917518a. --- yambopy/bse/realSpace_excitonwf.py | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 8bb59aa0..7731d44d 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -176,15 +176,6 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, assert nk == nkBZ, "kpoint mismatch" # fixed_postion = np.array(fixed_postion) - fixed_postion += np.array(supercell)//2 - # - if fix_particle == 'h': - print("Position of the hole is set to (reduced units) : ", - fixed_postion[0], fixed_postion[1], fixed_postion[2]) - elif fix_particle == 'e': - print("Position of the electron is set to (reduced units) : ", - fixed_postion[0], fixed_postion[1], fixed_postion[2]) - # # hole_bnds = [bse_bnds[0],bse_bnds[0]+nv] elec_bnds = [bse_bnds[0]+nv,bse_bnds[1]] @@ -223,7 +214,21 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # Compute nstates, nk, Nx, Ny, Nz object if out_res is None : print("Wfc FFT Grid : ",fft_box[0], fft_box[1], fft_box[2]) # - ktree = wfcdb.ktree + ## + # find the nearest fft grid point. + fx_pnt_int = np.floor(fixed_postion) + fixed_postion -= fx_pnt_int + fixed_postion = np.round(fixed_postion * fft_box) / fft_box + fixed_postion += fx_pnt_int + # shift the position of hole to middle of supercell + fixed_postion += np.array(supercell)//2 + # + if fix_particle == 'h': + print("Position of the hole is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) + if fix_particle == 'e': + print("Position of the electron is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) + # + ktree = wfcdb.ktree #build_ktree(wfcdb.kBZ) # nspinorr = wfcdb.nspinor if out_res is not None: From d9f39591bc4f8ae4dfd5f356329246f80e9c8bcc Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 14 Jul 2025 12:23:14 +0200 Subject: [PATCH 57/73] save progress --- yambopy/bse/realSpace_excitonwf.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 7731d44d..f7db539e 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -198,7 +198,7 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, for ik in range(len(wfcdb.gvecs)): idx_gvecs_tmp = np.arange(wfcdb.ngvecs[ik],dtype=int) if wfcCutoffRy > 0: - tmp_gvecs = np.linalg.norm((wfcdb.gvecs[ik, :wfcdb.ngvecs[ik], :] + tmp_gvecs = 2*np.pi*np.linalg.norm((wfcdb.gvecs[ik, :wfcdb.ngvecs[ik], :] + wfcdb.kpts_iBZ[ik][None,:])@blat,axis=-1) idx_tmp = tmp_gvecs < np.sqrt(wfcCutoffRy) idx_gvecs_tmp = idx_gvecs_tmp[idx_tmp].copy() @@ -224,9 +224,11 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, fixed_postion += np.array(supercell)//2 # if fix_particle == 'h': - print("Position of the hole is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) + print("Position of the hole (reduced units) is set to : ", + fixed_postion[0], fixed_postion[1], fixed_postion[2]) if fix_particle == 'e': - print("Position of the electron is set to : ", fixed_postion[0], fixed_postion[1], fixed_postion[2]) + print("Position of the electron (reduced units) is set to : ", + fixed_postion[0], fixed_postion[1], fixed_postion[2]) # ktree = wfcdb.ktree #build_ktree(wfcdb.kBZ) # From a6d7bad99ef6a79d986a3b9156a81d9d775613ed Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 21 Aug 2025 23:25:54 +0200 Subject: [PATCH 58/73] add spin support to all exciton stuff --- yambopy/bse/exciton_matrix_elements.py | 25 +++---- yambopy/bse/exciton_spin.py | 2 +- yambopy/bse/realSpace_excitonwf.py | 56 +++++++--------- yambopy/bse/rotate_excitonwf.py | 40 +++++------ yambopy/dbs/excitondb.py | 21 +++--- yambopy/symmetries/__init__.py | 14 ++++ yambopy/symmetries/crystal_symmetries.py | 84 ++++++++++++++++++++++++ 7 files changed, 168 insertions(+), 74 deletions(-) create mode 100644 yambopy/symmetries/__init__.py create mode 100644 yambopy/symmetries/crystal_symmetries.py diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index eed233c8..d58557df 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -26,12 +26,13 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di O_qvec : array_like Momentum transfer vector q in crystal coordinates. Akq : array_like - Wavefunction coefficients for k+q (bra wfc) with shape (n_exe_states, nk, nc, nv). + 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, nk, nc, nv). + 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, m_bnd, n_bnd). - ie Omn = < k+q, m | O(q) | n, k>, where m_bnd and n_bnd are final and initial bands, respectively. + 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 @@ -53,9 +54,9 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di # Number of arbitrary parameters (lambda) in the Omn matrix nlambda = Omn.shape[0] # - assert len(Akq.shape) == 4, "Works only with TDA." + assert Akq.shape[1] == 1, "Works only with TDA." # Shape of the wavefunction coefficients - n_exe_states, nk, nc, nv = Akq.shape + 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" @@ -64,19 +65,19 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di 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 not ktree : ktree = build_ktree(kpts) + 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:] # Occupied part - Ovv = Omn[:, idx_k_minus_Q_minus_q, :nv, :nv] # conduction part + 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[:, idx_k_minus_q, ...]) - Akq_conj = Akq.reshape(n_exe_states, -1).conj() + 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: @@ -88,7 +89,7 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di for il in range(nlambda): # Compute the electron contribution if contribution == 'e' or contribution == 'b': - tmp_wfc = Occ[il][None, :, :, :] @ Ak_electron + tmp_wfc = Occ[il][None, ...] @ Ak_electron # # Compute the hole contribution and subtract from the electron contribution if contribution == 'h' or contribution == 'b': diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 8f40188d..6785d3d2 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -62,7 +62,7 @@ def compute_exciton_spin(lattice, excdb, wfdb, elec_sz, contribution='b',diagona # # Compute the exciton spin matrix elements exe_Sz = exciton_X_matelem(excQpt, np.array([0, 0, 0]), Akcv, - Akcv, elec_sz[None, ...], wfdb.kBZ, + Akcv, elec_sz[None,:,None,...], wfdb.kBZ, diagonal_only=diagonal,contribution=contribution) # return exe_Sz[0] diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index f7db539e..70fa32b3 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -78,27 +78,19 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, - exe_wfc_real: Real-space exciton wavefunction (nstates, nspinor_electron, nspinor_hole, Nx_grid, Ny_grid, Nz_grid) """ - if len(Akcv.shape) == 5: - #nonTDA - ## first the resonat part - supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ + ## first the resonat part + supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ ex_wf2Real_kernel(Akcv[:,0], Qpt, wfcdb, bse_bnds, fixed_postion, - fix_particle=fix_particle, supercell=supercell, - wfcCutoffRy=wfcCutoffRy, block_size=block_size, + fix_particle=fix_particle, supercell=supercell, + wfcCutoffRy=wfcCutoffRy, block_size=block_size, ares=False, out_res=None) - # add to antiresonant part + # for nonTDA add the anti-resonant part + if Akcv.shape[1] == 2: supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ ex_wf2Real_kernel(Akcv[:,1], Qpt, wfcdb, bse_bnds, fixed_postion, fix_particle=fix_particle, supercell=supercell, wfcCutoffRy=wfcCutoffRy, block_size=block_size, ares=True, out_res=exe_wfc_real) - else: - assert len(Akcv.shape) == 4, "wrong Akcv dimensions" - supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ - ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, - fix_particle=fix_particle, supercell=supercell, - wfcCutoffRy=wfcCutoffRy, block_size=block_size, - ares=False, out_res=None) return supercell_latvecs,atom_nums,atom_pos,exe_wfc_real @@ -118,7 +110,7 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, Note: For density, one must compute the absolute value squared. Args: - Akcv (numpy.ndarray): Exciton coefficients [nstates,nk,nc,nv] + Akcv (numpy.ndarray): Exciton coefficients [nstates,ns,nk,nc,nv] Qpt (numpy.ndarray): Q-point in crystal coordinates [3] wfcdb (YamboWFDB): Wavefunction database bse_bnds (list): BSE band range used in bse [nb1, nb2]. fortran indexing. i.e index starts from 1. @@ -167,7 +159,8 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, "%d is used in bse but not found in wfcdb, load more wfcs" %(wfcdb.min_bnd + wfcdb.nbands) bse_bnds = np.array(bse_bnds,dtype=int)-wfcdb.min_bnd - nstates, nk, nc, nv = Akcv.shape + # ns is number of collinear spin components + nstates, ns, nk, nc, nv = Akcv.shape kpt_idx = wfcdb.ydb.kpoints_indexes sym_idx = wfcdb.ydb.symmetry_indexes @@ -183,11 +176,11 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, if ares: # if anti-resonant part, we swap the electron and hole bands # first transpose c,v dimensions - Akcv = Akcv.transpose(0,1,3,2) + Akcv = Akcv.transpose(0,1,2,4,3) tmp = hole_bnds hole_bnds = elec_bnds elec_bnds = tmp - nstates, nk, nc, nv = Akcv.shape + nstates, ns, nk, nc, nv = Akcv.shape lat_vec = wfcdb.ydb.lat.T blat = np.linalg.inv(lat_vec) @@ -234,12 +227,12 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # nspinorr = wfcdb.nspinor if out_res is not None: - exe_wfc_real = out_res.reshape(nstates, nspinorr, nspinorr, + exe_wfc_real = out_res.reshape(nstates, ns, nspinorr, nspinorr, supercell[0],fft_box[0], supercell[1],fft_box[1], supercell[2],fft_box[2]) else: - exe_wfc_real = np.zeros((nstates, nspinorr, nspinorr, + exe_wfc_real = np.zeros((nstates, ns, nspinorr, nspinorr, supercell[0],fft_box[0], supercell[1],fft_box[1], supercell[2],fft_box[2]), @@ -260,7 +253,8 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, FFFboxs = np.zeros((fft_box[0],fft_box[1],fft_box[2],3)) FFFboxs[...,0], FFFboxs[...,1], FFFboxs[...,2] = np.meshgrid(FFTxx,FFTyy,FFTzz,indexing='ij') # - exe_tmp_wf = np.zeros((nstates, nspinorr, nspinorr, min(nk,block_size) , + # ns is nspin + exe_tmp_wf = np.zeros((nstates, ns, nspinorr, nspinorr, min(nk,block_size) , fft_box[0], fft_box[1], fft_box[2]),dtype=np.complex64) # exp_tmp_kL = np.zeros((min(nk,block_size) ,supercell[0], @@ -343,11 +337,11 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, ft_wfcr *= exp_kx_r[None,None,None,...] # if fix_particle == 'h': - np.einsum('ncv,vy,cxijk->nxyijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], - optimize=True,out=exe_tmp_wf[:,:,:,ik-ikstart]) + np.einsum('nscv,svy,scxijk->nsxyijk',Akcv[:,:,ik,...],fx_wfc,ft_wfcr, + optimize=True,out=exe_tmp_wf[:,:,:,:,ik-ikstart]) else : - np.einsum('ncv,cx,vyijk->nxyijk',Akcv[:,ik,...],fx_wfc[0],ft_wfcr[0], - optimize=True,out=exe_tmp_wf[:,:,:,ik-ikstart]) + np.einsum('nscv,scx,svyijk->nsxyijk',Akcv[:,:,ik,...],fx_wfc,ft_wfcr, + optimize=True,out=exe_tmp_wf[:,:,:,:,ik-ikstart]) # #exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] exp_tmp_kL[ik-ikstart] = np.exp(1j*2*np.pi*np.einsum('...x,x->...',Lsupercells,ft_kvec)) @@ -356,19 +350,19 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, pbar.update(1) # ## perform gemm operation - total_gemms_t = nstates*nspinorr**2 + total_gemms_t = nstates*ns*nspinorr**2 exp_tmp_kL_tmp = exp_tmp_kL.reshape(len(exp_tmp_kL),-1)[:(ikstop-ikstart)].T - exe_tmp_wf_tmp = exe_tmp_wf.reshape(nstates,nspinorr,nspinorr,-1,np.prod(fft_box)) + exe_tmp_wf_tmp = exe_tmp_wf.reshape(nstates,ns,nspinorr,nspinorr,-1,np.prod(fft_box)) exe_tmp_wf_tmp = exe_tmp_wf_tmp[...,:(ikstop-ikstart),:] # for igemms in range(total_gemms_t): - ii, jj, kk = np.unravel_index(igemms, (nstates,nspinorr,nspinorr)) - Ctmp = (exp_tmp_kL_tmp @ exe_tmp_wf_tmp[ii, jj, kk ]) + ii, iis, jj, kk = np.unravel_index(igemms, (nstates,ns,nspinorr,nspinorr)) + Ctmp = (exp_tmp_kL_tmp @ exe_tmp_wf_tmp[ii, iis, jj, kk ]) Ctmp = Ctmp.reshape(supercell[0], supercell[1], supercell[2], fft_box[0], fft_box[1], fft_box[2]) - exe_wfc_real[ii, jj, kk ] += Ctmp.transpose(0,3,1,4,2,5) + exe_wfc_real[ii, iis, jj, kk ] += Ctmp.transpose(0,3,1,4,2,5) # - exe_wfc_real = exe_wfc_real.reshape(nstates,nspinorr,nspinorr, + exe_wfc_real = exe_wfc_real.reshape(nstates,ns,nspinorr,nspinorr, supercell[0]*fft_box[0], supercell[1]*fft_box[1], supercell[2]*fft_box[2]) diff --git a/yambopy/bse/rotate_excitonwf.py b/yambopy/bse/rotate_excitonwf.py index 7595aa97..0d548346 100644 --- a/yambopy/bse/rotate_excitonwf.py +++ b/yambopy/bse/rotate_excitonwf.py @@ -2,6 +2,7 @@ from yambopy.kpoints import build_ktree, find_kpt from yambopy.tools.function_profiler import func_profile + @func_profile def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=None): """ @@ -15,8 +16,8 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non Parameters ---------- Ak : array_like - Exciton wavefunction coefficients with shape (n_exe_states, nk, nc, nv). - In case of non-TDA, it is (n_exe_states, 2, nk, nc, nv) + Exciton wavefunction coefficients with shape (n_exe_states, 1 or 2, nspin, nk, nc, nv). + 1 for TDA and 2 for coupling symm_mat_red : array_like Symmetry matrix in reduced coordinates with shape (3, 3). kpoints : array_like @@ -24,7 +25,7 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non exe_qpt : array_like Momentum of the exciton (q-point) in crystal coordinates with shape (3,). dmats : array_like - Representation matrices for the symmetry operation with shape (nk, Rk_band, k_band). + Representation matrices for the symmetry operation with shape (nk, nspin, Rk_band, k_band). time_rev : bool If True, apply time-reversal symmetry to the wavefunction. ktree : object, optional @@ -39,12 +40,11 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non rot_Ak = np.zeros(Ak.shape, dtype=Ak.dtype) # Check TDA tda = True - if len(Ak.shape) != 4 : tda = False + if Ak.shape[1] == 2 : tda = False - nk, nc, nv = Ak.shape[-3:] + ns, nk, nc, nv = Ak.shape[2:] # Build a k-point tree if not provided - if not ktree: - ktree = build_ktree(kpoints) + if ktree is None: ktree = build_ktree(kpoints) # Compute the indices of Rk and Rk - q Rkpts = kpoints @ symm_mat_red.T # Rotated k-points @@ -53,24 +53,24 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non idx_k_minus_q = find_kpt(ktree, k_minus_q) # Indices of k - q # Extract the conduction and valence parts of the representation matrices - Dcc = dmats[:, nv:, nv:] # Conduction band part - Dvv = dmats[idx_k_minus_q, :nv, :nv].conj() # Valence band part (conjugated) + Dcc = dmats[:, :, nv:, nv:].transpose(1,0,2,3) # Conduction band part + Dvv = dmats[idx_k_minus_q, :, :nv, :nv].transpose(1,0,2,3).conj() # Valence band part (conjugated) # Apply time-reversal symmetry if required Ak_tmp = Ak - if time_rev: - Ak_tmp = Ak.conj() + if time_rev: Ak_tmp = Ak.conj() # Rotate the Ak wavefunction using the representation matrices - if tda : - rot_Ak[:, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp) @ (Dvv.transpose(0, 2, 1)[None, ...]) - else : - ## rotate the resonant part - rot_Ak[:, 0, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp[:,0,...]) @ (Dvv.transpose(0, 2, 1)[None, ...]) + ## rotate the resonant part + rot_Ak[:, :1, :, idx_Rk, ...] = ((Dcc[None, ...] @ Ak_tmp[:,0,...]) + @ (Dvv.transpose(0, 1, 3, 2)[None, ...]) + ).reshape(rot_Ak[:,:1].shape) + if not tda : # Rotate the anti-resonant part - Dvv = dmats[:, :nv, :nv] - Dcc = dmats[idx_k_minus_q, nv:, nv:].conj() - rot_Ak[:, 1, idx_Rk, ...] = (Dcc[None, ...] @ Ak_tmp[:,1,...]) @ (Dvv.transpose(0, 2, 1)[None, ...]) - # done return the rotate wfc + Dvv = dmats[:, :, :nv, :nv].transpose(1,0,2,3) + Dcc = dmats[idx_k_minus_q, :, nv:, nv:].transpose(1,0,2,3).conj() + rot_Ak[:, 1:, :, idx_Rk, ...] = ((Dcc[None, ...] @ Ak_tmp[:,1,...]) + @ (Dvv.transpose(0, 1, 3, 2)[None, ...]) + ).reshape(rot_Ak[:, 1:].shape) return rot_Ak diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index c31c3dc8..90a5fc4d 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -239,11 +239,11 @@ def write_sorted(self,prefix='yambo'): def get_Akcv(self): """ - Convert eigenvectors from (neigs,BS_table) -> (neigs,k,c,v) - For now, only works for nspin = 1. nspinor = 1/2 also works. + Convert eigenvectors from (neigs,BS_table) -> (neigs,nblks,nspin,k,c,v) + nblks = 2 for coupling, else 1 for TDA """ - - assert self.spin_pol == 'no', "Rearrange_Akcv works only for nspin = 1" + nspin = 1 + if self.spin_pol == 'pol': nspin = 2 # tmp_akcv = getattr(self, 'Akcv', None) if tmp_akcv is not None: return tmp_akcv @@ -254,19 +254,20 @@ def get_Akcv(self): nk = self.nkpoints nv = self.nvbands nc = self.ncbands - # Make sure nc * nv * nk = BS_TABLE length - table_len = nk*nv*nc - assert table_len == self.table.shape[0], "BS_TABLE length not equal to nc * nv * nk" + # Make sure nspin * nc * nv * nk = BS_TABLE length + table_len = nspin*nk*nv*nc + assert table_len == self.table.shape[0], "BS_TABLE length not equal to ns * nc * nv * nk" # v_min = np.min(self.table[:,1]) c_min = np.min(self.table[:,2]) bs_table0 = self.table[:,0]-1 bs_table1 = self.table[:,1] - v_min bs_table2 = self.table[:,2] - c_min + bs_table3 = self.table[:,3]-1 # eig_wfcs_returned = np.zeros(eig_wfcs.shape,dtype=eig_wfcs.dtype) # - sort_idx = bs_table0*nc*nv + bs_table2*nv + bs_table1 + sort_idx = bs_table0*nc*nv + bs_table2*nv + bs_table1 + nk*nc*nv*bs_table3 # eig_wfcs_returned[:,sort_idx] = eig_wfcs[...,:table_len] # check if this is coupling . @@ -274,9 +275,9 @@ def get_Akcv(self): eig_wfcs_returned[:,sort_idx+table_len] = eig_wfcs[...,table_len:] # NM : Note that here v and c are inverted i.e # psi_S = Akcv * phi_v(r_e) * phi_c^*(r_h) - eig_wfcs_returned = eig_wfcs_returned.reshape(-1,2,nk,nc,nv) + eig_wfcs_returned = eig_wfcs_returned.reshape(-1,2,nspin,nk,nc,nv) else : - eig_wfcs_returned = eig_wfcs_returned.reshape(-1,nk,nc,nv) + eig_wfcs_returned = eig_wfcs_returned.reshape(-1,1,nspin,nk,nc,nv) # self.Akcv = eig_wfcs_returned return self.Akcv diff --git a/yambopy/symmetries/__init__.py b/yambopy/symmetries/__init__.py new file mode 100644 index 00000000..5d8a8655 --- /dev/null +++ b/yambopy/symmetries/__init__.py @@ -0,0 +1,14 @@ +# Copyright (C) 2023, Claudio Attaccalite +# All rights reserved. +# +# This file is part of yambopy +# +""" +submodule with classes to handle post-processing observables related to spectroscopy +""" + +from .base_optical import BaseOpticalProperties +from .exciton_group_theory import ExcitonGroupTheory +from .ex_dipole import ExcitonDipole +from .ex_phonon import ExcitonPhonon +from .luminescence import Luminescence diff --git a/yambopy/symmetries/crystal_symmetries.py b/yambopy/symmetries/crystal_symmetries.py new file mode 100644 index 00000000..e7e6bb69 --- /dev/null +++ b/yambopy/symmetries/crystal_symmetries.py @@ -0,0 +1,84 @@ +import numpy as np +import spgrep +import spglib +from yambopy.dbs.latticedb import YamboLatticeDB + +class Crystal_Symmetries: + def __init__(self, latdb, magnetic_moments=None, tol=1e-5): + """ + Class to handle crystal symmetries + # + Given a lattice database, finds all symmetries of a crystal. + Works irrespective of whether the SAVE has symmetries or not. + # + Parameters: + ----------- + latdb : object + Lattice database object with required attributes: + magnetic_moments : array (natom) or (natom,3) or None, optional + magnetic_moments of each atom (default: None) + In collinear case, natom is sufficient, but in non-collinear + case provide (mx, my, mz) for each atom. + tol : float, optional + Symmetry tolerance (default: 1e-5) + Members: + -------- + dataset : dict + The full symmetry dataset from spglib + spacegroup_type : dict + Space group type information + rotations : ndarray + Rotation matrices in Cartesian coordinates (n_sym, 3, 3) + translations : ndarray + Translation vectors in Cartesian coordinates (n_sym, 3) + international_symbol : str + International symbol of the space group + hall_symbol : str + Hall symbol of the space group + wyckoffs : list + Wyckoff letters for each atom + pointgroup : str + Point group symbol + spacegroup_number : int + International space group number + """ + # Get lattice and atomic information + lattice = latdb.lat + numbers = np.rint(latdb.atomic_numbers).astype(int) + lat_inv = np.linalg.inv(lattice) + positions = latdb.car_atomic_positions @ lat_inv + positions = positions - np.floor(positions) + # + self.magnetic = not (magnetic_moments is None) + # + # Get symmetry dataset + if self.magnetic : + cell = (lattice, positions, numbers, magnetic_moments) + print("Currently Magnetic systems not supported.") + exit() + else : + cell = (lattice, positions, numbers) + self.dataset = spglib.get_symmetry_dataset(cell, symprec=tol) + self.spacegroup_type = spglib.get_spacegroup_type(self.dataset.hall_number) + # + # Convert rotation and fractional translations to Cartesian units + self.rotations = lattice.T[None,:,:] @ self.dataset.rotations @ lat_inv.T[None,:,:] + self.translations = self.dataset.translations @ lattice + # + self.international_symbol = self.dataset.international + self.hall_symbol = self.dataset.hall + self.wyckoffs = self.dataset.wyckoffs + self.pointgroup = self.dataset.pointgroup + self.pointgroup_schoenflies = self.spacegroup_type.pointgroup_schoenflies + self.spacegroup_number = self.dataset.number + + +## Test +if __name__ == "__main__": + import os + np.set_printoptions(suppress=True) + latdb = YamboLatticeDB.from_db_file(os.path.join('.', 'SAVE')) + symm = Crystal_Symmetries(latdb,tol=1e-4) + print(symm.wyckoffs) + print(symm.pointgroup_schoenflies) + From b8ede226bd0452139ce3fece948e0d1047f5cadb Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 21 Aug 2025 23:29:47 +0200 Subject: [PATCH 59/73] fix bug --- yambopy/dbs/excitondb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index 90a5fc4d..852ec5fd 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -337,7 +337,7 @@ def real_wf_to_cube(self, iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], density *= phase # # sum over spinor indices and degenerate states - real_wfc = np.sum(density,axis=(0,1,2)) + real_wfc = np.sum(density,axis=(0,1,2,3)) # normalize with max value max_normalize_val = np.max(np.abs(real_wfc)) print('Max Normalization value: ',max_normalize_val) From dbd59984a6bd8b9e0a03c6cf29642c3f4ee70c75 Mon Sep 17 00:00:00 2001 From: Murali Date: Thu, 21 Aug 2025 23:31:32 +0200 Subject: [PATCH 60/73] fix bug --- yambopy/symmetries/__init__.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/yambopy/symmetries/__init__.py b/yambopy/symmetries/__init__.py index 5d8a8655..22b7d0c7 100644 --- a/yambopy/symmetries/__init__.py +++ b/yambopy/symmetries/__init__.py @@ -1,14 +1,5 @@ -# Copyright (C) 2023, Claudio Attaccalite +# Copyright (C) 2025, YamboPy project # All rights reserved. # # This file is part of yambopy # -""" -submodule with classes to handle post-processing observables related to spectroscopy -""" - -from .base_optical import BaseOpticalProperties -from .exciton_group_theory import ExcitonGroupTheory -from .ex_dipole import ExcitonDipole -from .ex_phonon import ExcitonPhonon -from .luminescence import Luminescence From 92e418ee26cfec1b3c3f47c59df9f81372d8aac8 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 22 Aug 2025 10:06:20 +0200 Subject: [PATCH 61/73] remove comment --- yambopy/bse/realSpace_excitonwf.py | 1 - 1 file changed, 1 deletion(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 70fa32b3..82e21e6d 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -328,7 +328,6 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, exp_fx = np.exp(2*np.pi*1j*((fx_gvec + fx_kvec[None,:])@fixed_postion)) fx_wfc *= exp_fx[None,None,None,:] fx_wfc = np.sum(fx_wfc,axis=-1) #(spin,bnd,spinor) - ## for now only nspin = 1 works. ns1, nbndc, nspinorr, ng = ft_wfc.shape #if ft_ikpt not in prev_ikpts: ft_wfcr = wfcdb.to_real_space(ft_wfc.reshape(-1,nspinorr,ng),ft_gvec, grid=fft_box) From a1e764ba964aa532a04bf9ca1d9e4e5cf87f5e24 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 22 Aug 2025 10:09:11 +0200 Subject: [PATCH 62/73] add comment --- yambopy/bse/realSpace_excitonwf.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 82e21e6d..62bed7c9 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -75,8 +75,9 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, - supercell_latvecs: Supercell lattice vectors (3,3) - atom_nums: Atomic numbers. - atom_pos: Atomic positions in cartisian units - - exe_wfc_real: Real-space exciton wavefunction (nstates, nspinor_electron, nspinor_hole, + - exe_wfc_real: Real-space exciton wavefunction (nstates, nspin, nspinor_electron, nspinor_hole, Nx_grid, Ny_grid, Nz_grid) + note if nspin =2, then nspinor = 1, similarly, if nspinor = 2, nspin = 1. """ ## first the resonat part supercell_latvecs,atom_nums,atom_pos,exe_wfc_real = \ @@ -135,7 +136,8 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, - atom_nums (numpy.ndarray): Atomic numbers [Natoms] - atom_pos (numpy.ndarray): Atomic positions in cartisian units [Natoms,3] - exe_wfc_real (numpy.ndarray): Wavefunction in real space - [nstates, nspinor_electron, nspinor_hole, FFTx, FFTy, FFTz] + [nstates, nspin, nspinor_electron, nspinor_hole, FFTx, FFTy, FFTz] + note if nspin =2, then nspinor = 1, similarly, if nspinor = 2, nspin = 1 """ # # From 401683aedc2fffdba6910ed34f68a616035f6072 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 22 Aug 2025 10:16:07 +0200 Subject: [PATCH 63/73] add comment --- yambopy/bse/realSpace_excitonwf.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 62bed7c9..901a0878 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -58,8 +58,8 @@ def ex_wf2Real(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, Args: Akcv (numpy.ndarray): Exciton wavefunction coefficients with shape: - - (Nstates,k,c,v) for TDA - - (Nstates,2,k,c,v) for non-TDA (2 for resonant/anti-resonant) + - (Nstates,1,ns,k,c,v) for TDA + - (Nstates,2,ns,k,c,v) for non-TDA (2 for resonant/anti-resonant) Qpt (numpy.ndarray): Q-point of exciton in crystal coordinates wfcdb (YamboWFDB): Wavefunction database bse_bnds (list): Band range used in BSE [min_band, max_band] (python indexing) From 0b4e10ba18e0655801d78186cea46d6a18b1672e Mon Sep 17 00:00:00 2001 From: Murali Date: Sat, 30 Aug 2025 08:21:05 +0200 Subject: [PATCH 64/73] steven's fix --- yambopy/bse/realSpace_excitonwf.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 901a0878..87d41a40 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -211,6 +211,7 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, # ## # find the nearest fft grid point. + fixed_postion = fixed_postion.astype(np.float64) fx_pnt_int = np.floor(fixed_postion) fixed_postion -= fx_pnt_int fixed_postion = np.round(fixed_postion * fft_box) / fft_box @@ -337,11 +338,14 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, exp_kx_r = np.exp(2*np.pi*1j*FFFboxs.reshape(-1,3)@ft_kvec).reshape(FFFboxs.shape[:3]) ft_wfcr *= exp_kx_r[None,None,None,...] # + fx_wfc = fx_wfc.astype(np.complex64) + ft_wfcr = ft_wfcr.astype(np.complex64) + # if fix_particle == 'h': - np.einsum('nscv,svy,scxijk->nsxyijk',Akcv[:,:,ik,...],fx_wfc,ft_wfcr, + np.einsum('nscv,svy,scxijk->nsxyijk',Akcv[:,:,ik,...].astype(np.complex64),fx_wfc,ft_wfcr, optimize=True,out=exe_tmp_wf[:,:,:,:,ik-ikstart]) else : - np.einsum('nscv,scx,svyijk->nsxyijk',Akcv[:,:,ik,...],fx_wfc,ft_wfcr, + np.einsum('nscv,scx,svyijk->nsxyijk',Akcv[:,:,ik,...].astype(np.complex64),fx_wfc,ft_wfcr, optimize=True,out=exe_tmp_wf[:,:,:,:,ik-ikstart]) # #exe_tmp_wf[:,:,:,ik-ikstart] *= exp_kx_r[...].reshape(FFFboxs.shape[:3])[None,None,None] From 77d22384be3b3b13e3519478422ce0c220bbe683 Mon Sep 17 00:00:00 2001 From: Murali Date: Wed, 10 Sep 2025 15:14:45 +0200 Subject: [PATCH 65/73] overlap wannier --- yambopy/dbs/wfdb.py | 69 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 50fdea03..7cc27f06 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -617,6 +617,75 @@ def Dmat(self, symm_mat=None, frac_vec=None, time_rev=None): # return Dmat + def OverlapUkkp(self, kpt_bra, kpt_ket): + """ + Compute the following matrix elements : < k_bra | e^{i(k_bra-k_ket).r} | k_ket> + in other words, it computes overlap of periodic parts of k_bra an k_ket + """ + kpt_bra = np.array(kpt_bra) + kpt_ket = np.array(kpt_ket) + + ikpt_ket = find_kpt(self.ktree,kpt_ket) + ikpt_bra = find_kpt(self.ktree,kpt_bra) + # + kpt_idx = self.ydb.kpoints_indexes + sym_idx = self.ydb.symmetry_indexes + # + ibz_ket = kpt_idx[ikpt_ket] + isym_ket = sym_idx[ikpt_ket] + # + ibz_bra = kpt_idx[ikpt_bra] + isym_bra = sym_idx[ikpt_bra] + # + ## get the wfcs: + k_rk_ket, w_rk_ket, g_rk_ket = self.rotate_wfc(ibz_ket, isym_ket) + k_rk_bra, w_rk_bra, g_rk_bra = self.rotate_wfc(ibz_bra, isym_bra) + + G0_ket = kpt_ket-k_rk_ket + G0_bra = kpt_bra-k_rk_bra + # + G0 = G0_ket-G0_bra + return wfc_inner_product(G0, w_rk_bra, g_rk_bra, np.array([0,0,0]), w_rk_ket, g_rk_ket) + + #def OverlapUkkp(self, kpt_bra, kpt_ket): + # """ + # Compute the following matrix elements : < k_bra | e^{i(k_bra-k_ket).r} | k_ket> + # in other words, it computes overlap of periodic parts of k_bra an k_ket + # """ + # kpt_bra = np.array(kpt_bra) + # kpt_ket = np.array(kpt_ket) + + # ikpt_ket = find_kpt(self.ktree,kpt_ket) + # ikpt_bra = find_kpt(self.ktree,kpt_bra) + # # + # kpt_idx = self.ydb.kpoints_indexes + # sym_idx = self.ydb.symmetry_indexes + # # + # ibz_ket = kpt_idx[ikpt_ket] + # isym_ket = sym_idx[ikpt_ket] + # # + # ibz_bra = kpt_idx[ikpt_bra] + # isym_bra = sym_idx[ikpt_bra] + # # + # ## get the wfcs: + # k_rk_ket, w_rk_ket, g_rk_ket = self.rotate_wfc(ibz_ket, isym_ket) + # k_rk_bra, w_rk_bra, g_rk_bra = self.rotate_wfc(ibz_bra, isym_bra) + + # grid = [30,30,180] + # wf_b = self.to_real_space(w_rk_bra[0], g_rk_bra, grid=grid) + # wf_k = self.to_real_space(w_rk_ket[0], g_rk_ket, grid=grid) + + # cel_vol = abs(np.linalg.det(self.ydb.lat.T)) + # fft_r = np.meshgrid(np.arange(grid[0])/grid[0],np.arange(grid[1])/grid[1],np.arange(grid[2])/grid[2],indexing='ij') + # rgrid = np.zeros((grid[0],grid[1],grid[2],3)) + # rgrid[...,0],rgrid[...,1],rgrid[...,2] = fft_r + # kdotr = np.einsum('...x,x->...',rgrid,(kpt_bra-kpt_ket)) + # return np.einsum('m...,n...->mn',np.conj(wf_b),wf_k*np.exp(1j*2*np.pi*kdotr),optimize=True)*cel_vol/np.prod(grid) + + + + + @func_profile def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gtree=None): """ From be0041344b06a1631b4ed6f3c6a3e2d5078c250e Mon Sep 17 00:00:00 2001 From: Murali Date: Wed, 10 Sep 2025 17:37:49 +0200 Subject: [PATCH 66/73] cleaup --- yambopy/dbs/wfdb.py | 40 ++-------------------------------------- 1 file changed, 2 insertions(+), 38 deletions(-) diff --git a/yambopy/dbs/wfdb.py b/yambopy/dbs/wfdb.py index 7cc27f06..93034d0b 100644 --- a/yambopy/dbs/wfdb.py +++ b/yambopy/dbs/wfdb.py @@ -647,45 +647,9 @@ def OverlapUkkp(self, kpt_bra, kpt_ket): G0 = G0_ket-G0_bra return wfc_inner_product(G0, w_rk_bra, g_rk_bra, np.array([0,0,0]), w_rk_ket, g_rk_ket) - #def OverlapUkkp(self, kpt_bra, kpt_ket): - # """ - # Compute the following matrix elements : < k_bra | e^{i(k_bra-k_ket).r} | k_ket> - # in other words, it computes overlap of periodic parts of k_bra an k_ket - # """ - # kpt_bra = np.array(kpt_bra) - # kpt_ket = np.array(kpt_ket) - - # ikpt_ket = find_kpt(self.ktree,kpt_ket) - # ikpt_bra = find_kpt(self.ktree,kpt_bra) - # # - # kpt_idx = self.ydb.kpoints_indexes - # sym_idx = self.ydb.symmetry_indexes - # # - # ibz_ket = kpt_idx[ikpt_ket] - # isym_ket = sym_idx[ikpt_ket] - # # - # ibz_bra = kpt_idx[ikpt_bra] - # isym_bra = sym_idx[ikpt_bra] - # # - # ## get the wfcs: - # k_rk_ket, w_rk_ket, g_rk_ket = self.rotate_wfc(ibz_ket, isym_ket) - # k_rk_bra, w_rk_bra, g_rk_bra = self.rotate_wfc(ibz_bra, isym_bra) - - # grid = [30,30,180] - # wf_b = self.to_real_space(w_rk_bra[0], g_rk_bra, grid=grid) - # wf_k = self.to_real_space(w_rk_ket[0], g_rk_ket, grid=grid) - - # cel_vol = abs(np.linalg.det(self.ydb.lat.T)) - # fft_r = np.meshgrid(np.arange(grid[0])/grid[0],np.arange(grid[1])/grid[1],np.arange(grid[2])/grid[2],indexing='ij') - # rgrid = np.zeros((grid[0],grid[1],grid[2],3)) - # rgrid[...,0],rgrid[...,1],rgrid[...,2] = fft_r - # kdotr = np.einsum('...x,x->...',rgrid,(kpt_bra-kpt_ket)) - # return np.einsum('m...,n...->mn',np.conj(wf_b),wf_k*np.exp(1j*2*np.pi*kdotr),optimize=True)*cel_vol/np.prod(grid) - - - - +## end of class +## @func_profile def wfc_inner_product(k_bra, wfc_bra, gvec_bra, k_ket, wfc_ket, gvec_ket, ket_Gtree=None): """ From 3414dcaa3631efa124ab2247fea98e6eb4445bdc Mon Sep 17 00:00:00 2001 From: Murali Date: Wed, 10 Sep 2025 19:53:41 +0200 Subject: [PATCH 67/73] fix wrong normalization for nonTDA --- yambopy/bse/realSpace_excitonwf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/bse/realSpace_excitonwf.py b/yambopy/bse/realSpace_excitonwf.py index 87d41a40..19821322 100644 --- a/yambopy/bse/realSpace_excitonwf.py +++ b/yambopy/bse/realSpace_excitonwf.py @@ -365,13 +365,13 @@ def ex_wf2Real_kernel(Akcv, Qpt, wfcdb, bse_bnds, fixed_postion, Ctmp = (exp_tmp_kL_tmp @ exe_tmp_wf_tmp[ii, iis, jj, kk ]) Ctmp = Ctmp.reshape(supercell[0], supercell[1], supercell[2], fft_box[0], fft_box[1], fft_box[2]) + Ctmp *= (1.0/np.prod(supercell)) exe_wfc_real[ii, iis, jj, kk ] += Ctmp.transpose(0,3,1,4,2,5) # exe_wfc_real = exe_wfc_real.reshape(nstates,ns,nspinorr,nspinorr, supercell[0]*fft_box[0], supercell[1]*fft_box[1], supercell[2]*fft_box[2]) - exe_wfc_real *= (1.0/np.prod(supercell)) # # compute postioon of fixed particle in cart units fixed_postion_cc = lat_vec@fixed_postion From 29eca3687864632811e7823c74abaecd05418292 Mon Sep 17 00:00:00 2001 From: Murali Date: Fri, 12 Sep 2025 15:32:47 +0200 Subject: [PATCH 68/73] fix bug --- yambopy/dbs/excitondb.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yambopy/dbs/excitondb.py b/yambopy/dbs/excitondb.py index 852ec5fd..b9df8c52 100644 --- a/yambopy/dbs/excitondb.py +++ b/yambopy/dbs/excitondb.py @@ -317,7 +317,7 @@ def real_wf_to_cube(self, iexe, wfdb, fixed_postion=[0,0,0], supercell=[1,1,1], if fix_particle == 'h': name_file = 'electron' else: name_file = 'hole' - if phase and real_wfc.shape[1] != 1: + if phase and wfdb.wf.shape[1] != 1 and wfdb.wf.shape[3] != 1: print("phase plot only works for nspin = 1 and nspinor == 1") phase = False if phase and len(iexe_degen_states) > 1: From 9ddd7b8531dc6c668f1c2f6e1b419e489a6f0bc2 Mon Sep 17 00:00:00 2001 From: Murali Date: Mon, 15 Sep 2025 23:02:23 +0200 Subject: [PATCH 69/73] spin support --- yambopy/bse/exciton_matrix_elements.py | 2 +- yambopy/letzelphc_interface/lelphcdb.py | 26 +++++++++++++++++++++---- 2 files changed, 23 insertions(+), 5 deletions(-) diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index d58557df..0a88c504 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -93,7 +93,7 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di # # Compute the hole contribution and subtract from the electron contribution if contribution == 'h' or contribution == 'b': - tmp_h = -Ak @ Ovv[il][None, ...] + tmp_h = -Ak[:,0] @ Ovv[il][None, ...] if contribution == 'b': tmp_wfc += tmp_h else: diff --git a/yambopy/letzelphc_interface/lelphcdb.py b/yambopy/letzelphc_interface/lelphcdb.py index 92f9d34f..d53d8ece 100644 --- a/yambopy/letzelphc_interface/lelphcdb.py +++ b/yambopy/letzelphc_interface/lelphcdb.py @@ -48,7 +48,9 @@ def __init__(self,filename,read_all=True,div_by_energies=True,verbose=False): self.nq = database.dimensions['nq'].size self.ns = database.dimensions['nspin'].size self.nsym = database.dimensions['nsym_ph'].size - + self.div_by_energies = div_by_energies # if true, the elph store are normalized with 1/(2*w_ph) + # + # conv = database['convention'][...].data if isinstance(conv, np.ndarray): if conv.dtype.kind == 'S': # Byte strings (C chars) @@ -59,12 +61,14 @@ def __init__(self,filename,read_all=True,div_by_energies=True,verbose=False): conv = conv.decode('utf-8').strip() else: conv = str(conv).strip() - stard_conv = False - if conv.strip() == 'standard': + conv = conv.strip().replace('\0', '') + # + # + if conv == 'standard': print("Convention used in Letzelphc : k -> k+q (standard)") else: print("Convention used in Letzelphc : k-q -> k (yambo)") - self.convention = conv.strip() + self.convention = conv # # Read DB self.kpoints = database.variables['kpoints'][:] @@ -182,6 +186,7 @@ def read_iq(self,iq, bands_range=[], database=None, convention='yambo'): The phonon eigenvectors. - ph_elph_me : ndarray The electron-phonon matrix elements with the specified convention. + ( nk, nm, nspin, initial bnd, final bnd) """ # if len(bands_range) == 0: @@ -209,7 +214,20 @@ def read_iq(self,iq, bands_range=[], database=None, convention='yambo'): ph_eigs = database['POLARIZATION_VECTORS'][iq,...].data eph_mat = eph_mat[...,0] + 1j*eph_mat[...,1] ph_eigs = ph_eigs[...,0] + 1j*ph_eigs[...,1] + ## normalize with ph_freq + if self.div_by_energies: + ph_freq_iq = np.sqrt(2.0*np.abs(self.ph_energies[iq])/(ha2ev/2.)) + if iq >0: + ph_freq_iq = 1.0/ph_freq_iq + eph_mat *= ph_freq_iq[None,:,None,None,None] + else: + eph_mat[:,:3] = 0.0 + ph_freq_iq = 1.0/ph_freq_iq[3:] + eph_mat[:,3:] *= ph_freq_iq[None,:,None,None,None] + if close_file :database.close() + ## output elph matrix elements unit (Ry if div_by_energies else Ry^1.5) + # ( nk, nm, nspin, initial bnd, final bnd) return [ph_eigs, self.change_convention(self.qpoints[iq],eph_mat, convention)] def change_convention(self, qpt, elph_iq, convention='yambo'): From e957ce77964baf6844c77879f8e05626c1d59a2b Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 16 Sep 2025 00:18:16 +0200 Subject: [PATCH 70/73] cleanup comments --- yambopy/bse/exciton_spin.py | 39 ++------------------------ yambopy/tools/degeneracy_finder.py | 45 +++++++++++++++--------------- 2 files changed, 25 insertions(+), 59 deletions(-) diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 6785d3d2..24cf757a 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -76,43 +76,8 @@ def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, sz=0.5 * np.array([[1, 0], [0, -1]]), return_dbs_and_spin=True): """ - Compute the spin matrix elements ⟨S'|S_z|S⟩ for excitons. - - This function calculates the spin matrix elements for excitons using - wavefunctions and spin operators. Easy to use interface. Use - compute_exciton_spin() incase you already have those db's - - Parameters - ---------- - path : str, optional - Path to the directory containing the `SAVE` folder. Default is `.`. - bse_dir : str, optional - Directory containing the BSE (Bethe-Salpeter Equation) data. Default is `'SAVE'`. - iqpt : int or list/array of ints, optional - Index or indices of the q-point(s) for which the exciton spin is computed. - Default is `1` (Gamma point). - nstates : int, optional - Number of exciton states to consider. If `-1`, all states are included. Default is `-1`. - contribution : {'b', 'e', 'h'}, optional - Specifies which contribution to compute: - - `'b'`: Total exciton spin (default). - - `'e'`: Electron spin only. - - `'h'`: Hole spin only. - degen_tol : float, optional - Degeneracy tolerance for excitons in eV. Default is `1e-2` eV. - sz : ndarray, optional - Spin-z operator matrix in the basis of spinor wavefunctions. - Default is `0.5 * np.array([[1, 0], [0, -1]])`. - return_dbs_and_spin : bool, optional - return [latticedb, wfdb, excdb (s), elec_spin_matrix] - Default is True - Returns - ------- - exe_Sz : ndarray - Spin matrix elements for excitons with shape `(niq, nstates, nstates)`, - where `niq` is the number of q-points. - if return_dbs_and_spin is True, will also return - [latticedb, wfdb, excdb (s), elec_spin_matrix] + Compute expectation value of S_z operator for excitons. + Examples -------- Compute the total spin matrix elements for excitons: diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py index f2b42f07..53f52788 100644 --- a/yambopy/tools/degeneracy_finder.py +++ b/yambopy/tools/degeneracy_finder.py @@ -48,32 +48,33 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): # NM : This is to ensure that in case the sequecence if # Arthematic progession with d < tol, we make sure that the # values are within its mean - final_degen_sets = [] - for group in degen_sets_sorted: - if len(group) == 0: - continue - group_eigenvalues = eigenvalues[group] - current_group = [group[0]] - current_mean = group_eigenvalues[0] + # final_degen_sets = [] + # for group in degen_sets_sorted: + # if len(group) == 0: + # continue + # group_eigenvalues = eigenvalues[group] + # current_group = [group[0]] + # current_mean = group_eigenvalues[0] - for i in range(1, len(group)): - diff = np.abs(group_eigenvalues[i] - current_mean) - tolerance = atol + rtol * np.abs(current_mean) + # for i in range(1, len(group)): + # diff = np.abs(group_eigenvalues[i] - current_mean) + # tolerance = atol + rtol * np.abs(current_mean) - if diff <= tolerance: - current_group.append(group[i]) - # Update the mean of the current group incrementally - current_mean = (current_mean * (len(current_group) - 1) + group_eigenvalues[i]) / len(current_group) - else: - final_degen_sets.append(np.array(current_group)) - current_group = [group[i]] - current_mean = group_eigenvalues[i] + # if diff <= tolerance: + # current_group.append(group[i]) + # # Update the mean of the current group incrementally + # current_mean = (current_mean * (len(current_group) - 1) + group_eigenvalues[i]) / len(current_group) + # else: + # final_degen_sets.append(np.array(current_group)) + # current_group = [group[i]] + # current_mean = group_eigenvalues[i] - # Append the last group - if current_group: - final_degen_sets.append(np.array(current_group)) + # # Append the last group + # if current_group: + # final_degen_sets.append(np.array(current_group)) - return final_degen_sets + # return final_degen_sets + return degen_sets_sorted if __name__ == '__main__': From 8278059658380b07a751d700351359962a032cb8 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 16 Sep 2025 10:23:02 +0200 Subject: [PATCH 71/73] update docs --- yambopy/bse/exciton_matrix_elements.py | 4 +-- yambopy/bse/exciton_spin.py | 36 +++++++++++++++++++++++++- yambopy/bse/rotate_excitonwf.py | 2 +- 3 files changed, 38 insertions(+), 4 deletions(-) diff --git a/yambopy/bse/exciton_matrix_elements.py b/yambopy/bse/exciton_matrix_elements.py index 0a88c504..65797a70 100644 --- a/yambopy/bse/exciton_matrix_elements.py +++ b/yambopy/bse/exciton_matrix_elements.py @@ -22,9 +22,9 @@ def exciton_X_matelem(exe_kvec, O_qvec, Akq, Ak, Omn, kpts, contribution='b', di Parameters ---------- exe_kvec : array_like - Exciton k-vector in crystal coordinates. + Exciton k-vector in crystal coordinates (k). O_qvec : array_like - Momentum transfer vector q in crystal coordinates. + 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 diff --git a/yambopy/bse/exciton_spin.py b/yambopy/bse/exciton_spin.py index 24cf757a..8509cc73 100644 --- a/yambopy/bse/exciton_spin.py +++ b/yambopy/bse/exciton_spin.py @@ -76,8 +76,42 @@ def compute_exc_spin_iqpt(path='.', bse_dir='SAVE', iqpt=1, sz=0.5 * np.array([[1, 0], [0, -1]]), return_dbs_and_spin=True): """ - Compute expectation value of S_z operator for excitons. + + 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: diff --git a/yambopy/bse/rotate_excitonwf.py b/yambopy/bse/rotate_excitonwf.py index 0d548346..8cac01ef 100644 --- a/yambopy/bse/rotate_excitonwf.py +++ b/yambopy/bse/rotate_excitonwf.py @@ -23,7 +23,7 @@ def rotate_exc_wf(Ak, symm_mat_red, kpoints, exe_qpt, dmats, time_rev, ktree=Non kpoints : array_like K-points in the full Brillouin zone (crystal coordinates) with shape (nk, 3). exe_qpt : array_like - Momentum of the exciton (q-point) in crystal coordinates with shape (3,). + Momentum of the given exciton (q-point) in crystal coordinates with shape (3,). dmats : array_like Representation matrices for the symmetry operation with shape (nk, nspin, Rk_band, k_band). time_rev : bool From 1b9806921b4793a93b0549150f38015a210be9a5 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 16 Sep 2025 10:27:09 +0200 Subject: [PATCH 72/73] update degenecy finder --- yambopy/tools/degeneracy_finder.py | 32 +----------------------------- 1 file changed, 1 insertion(+), 31 deletions(-) diff --git a/yambopy/tools/degeneracy_finder.py b/yambopy/tools/degeneracy_finder.py index 53f52788..d4e8d48d 100644 --- a/yambopy/tools/degeneracy_finder.py +++ b/yambopy/tools/degeneracy_finder.py @@ -43,37 +43,7 @@ def find_degeneracy_evs(eigenvalues, atol=1e-3, rtol=1e-3): # Group indices of degenerate states (in sorted order) degen_sets_sorted = np.split(idx_sorted, split_indices + 1) - - # Further split groups based on the mean of the group - # NM : This is to ensure that in case the sequecence if - # Arthematic progession with d < tol, we make sure that the - # values are within its mean - # final_degen_sets = [] - # for group in degen_sets_sorted: - # if len(group) == 0: - # continue - # group_eigenvalues = eigenvalues[group] - # current_group = [group[0]] - # current_mean = group_eigenvalues[0] - - # for i in range(1, len(group)): - # diff = np.abs(group_eigenvalues[i] - current_mean) - # tolerance = atol + rtol * np.abs(current_mean) - - # if diff <= tolerance: - # current_group.append(group[i]) - # # Update the mean of the current group incrementally - # current_mean = (current_mean * (len(current_group) - 1) + group_eigenvalues[i]) / len(current_group) - # else: - # final_degen_sets.append(np.array(current_group)) - # current_group = [group[i]] - # current_mean = group_eigenvalues[i] - - # # Append the last group - # if current_group: - # final_degen_sets.append(np.array(current_group)) - - # return final_degen_sets + return degen_sets_sorted From ddee1da72f8075ae739581369a654f2181ccb808 Mon Sep 17 00:00:00 2001 From: Murali Date: Tue, 16 Sep 2025 10:35:49 +0200 Subject: [PATCH 73/73] update tutorial for ex-wf --- tutorial/databases_yambopy/exc_real_wf.py | 40 +++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 tutorial/databases_yambopy/exc_real_wf.py diff --git a/tutorial/databases_yambopy/exc_real_wf.py b/tutorial/databases_yambopy/exc_real_wf.py new file mode 100644 index 00000000..aab96d62 --- /dev/null +++ b/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 ! +