Chirality in monolayer MoS2

[rishabhsrsw7]

Dear developers,

I have been trying to recreate the reciprocal space plot of monolayer MoSe_2 in https://arxiv.org/pdf/2511.21540 (Fig. 5(b)), for monolayer MoS_2. I have tried changing the code of the tutorial using the existing exciton k space plot https://github.com/yambo-code/yambopy/b ... exciton.py,

Code: Select all

#Authors: MN 
"""
Script to compute total crystal angular momentum
This script loads Yambo databases, calculates the angular momentum for a specific exciton,
and generates a reciprocal-space (K-space) visualization.
"""

import numpy as np
import os
import matplotlib.pyplot as plt
from yambopy.dbs.excitondb import YamboExcitonDB
from yambopy.dbs.latticedb import YamboLatticeDB
from yambopy.dbs.wfdb import YamboWFDB

## Inputs
iqpt = 1
path = '.'
gw_bse_dir = 'GW_BSE'
iexe = 2
degen_tol = 1e-4

# ======= Load data ========

# Load the lattice db
lattice = YamboLatticeDB.from_db_file(os.path.join(path, 'SAVE', 'ns.db1'))

# Load exciton db
filename = 'ndb.BS_diago_Q%d' % (iqpt)
excdb = YamboExcitonDB.from_db_file(lattice, filename=filename,
                                    folder=os.path.join(path, gw_bse_dir), neigs = iexe + 101)

# Load the wavefunction database
wfdb = YamboWFDB(path=path, latdb=lattice, bands_range=[np.min(excdb.table[:, 1]) - 1, np.max(excdb.table[:, 2])])

# Symmetry definitions
symm_mat_cart = lattice.sym_car[3]
frac_trans_cart = np.zeros(3)

# Compute Total Crystal Angular Momentum
sbasis = excdb.total_crys_angular_momentum(wfdb, iexe, symm_mat_cart, frac_trans_cart, degen_tol=degen_tol)

jvals = sbasis[0]
Akcv_j = sbasis[1]

print('Total crystal angular momentum : ', jvals)

# Find degenerate states (0-based indexing for internal arrays)
idegen = np.array(excdb.get_degenerate(iexe + 1, eps=degen_tol), dtype=int) - 1

# Update the existing wfs with the new rotated simultaneous eigenbasis
old_eig_states = excdb.Akcv[idegen].copy()
excdb.Akcv[idegen] = Akcv_j

# ======= Plotting in K-Space ========

print("="*80)
for i in range(len(jvals)):
    state_idx = idegen[i] + 1
    j_val = jvals[i]

    print('*********** Plotting K-space wf for state %d with j = %.2f ***************' % (state_idx, j_val))
    print("="*80)

    fig = plt.figure(figsize=(6,6))
    ax  = fig.add_axes([0.15, 0.15, 0.80, 0.80])

    excdb.plot_exciton_2D_ax(ax, [state_idx], mode='hexagon', limfactor=0.8, scale=180)

    ax.set_title(f"Exciton State {state_idx} | j = {j_val:.2f}")

    plt.savefig(f"Kspace_state_{state_idx}_j_{j_val:.2f}.png", dpi=300, bbox_inches='tight')
    plt.show()

print("="*80)

# Restore the original database coefficients to leave everything as it was
excdb.Akcv[idegen] = old_eig_states

.

The results are however disappointing and I am unable to understand if my code is wrong or there is a mistake I am making in the calculation. I am attaching the figures that this code generates, and if anybody can help me with this, it'd be greatly appreciated.

[palful]

Dear Rishabh,

You actually didn't do anything wrong, the problem is that the function `plot_exciton_2D_ax` that you call to plot in k-space, actually recalculates the weights using the unrotated basis, i.e., it calls the function `get_exciton_weights` that returns the |A_kcv|^2 values calculated from `self.eigenvectors` rather than `self.Akcv` which is the one that was overwritten after the angular momentum basis rotation.

The reason is simply historical since that bit of code was written earlier than the part about total crystal angular momentum which is a new implementation. In the future, it probably would be better if the function `plot_exciton_2D_ax` first checks for the existence of self.Akcv and uses it if found. I will note it for a future patch.

What you can do is simply to make your own scatter plot with matplotlib instead of relying on `plot_exciton_2D_ax`. You can basically copy the plotting part for that function, and calculate the weights |\sum_cv A_cv(k)|^2 yourself starting from excdb.Akcv. This should provide you with the results you want!

Best,
Fulvio

PS: you can also change the size of the markers in the scatterplot according to the k-mesh you are using, so as to avoid them overlapping.

[rishabhsrsw7]

Dear Fulvio,

Thanks a lot for your help with this, I think I was able to update the code correctly, and run it for monolayer hBN, where I can see the segregation for either circular polarization.

However, I wanted to ask if this is also affected by the databases we use for generating the rotated |A_kcv|^2. I tried the same for monolayer MoS_2 where I used SLEPC database, and for the first bright degenerate exciton, my weights are localised at the same position irrespective of the polarization. .

Regards
Rishabh

[palful]

Dear Rishabh,

It should work also in the case of the slepc solver, as long as:

- you correctly include the degenerate subspace of the "chiral" states, i.e., in your case of the bright A exciton it should be `iexe=2` (python index), then the code should find state `iexe=3` automatically.
- you correctly provide the C_3 rotation symmetry matrix `symm_mat_cart`, although since monolayer MoS2 and monolayer hBN have the same point group (if I am not mistaken) it should have the same index as in the previous case you did.

If you check these two things and still it does not work as intended, you may attach your script and I will try to test it over a small MoS2 calculation that I have.

Best,
Fulvio