shankar1729 / jdftx

JDFTx: software for joint density functional theory
http://jdftx.org
79 stars 49 forks source link

Eliashberg spectral function α2F in terms of the mode-resolved coupling strengths λqν #312

Closed dijasila closed 1 month ago

dijasila commented 5 months ago

Hello everyone, Is it possible in the current implementation to calculate the The Eliashberg spectral function α2F d in terms of the mode-resolved coupling strengths λqν and the phonon frequencies? is there a simple way of doing that? All I want is to precisely know which peak in the α2F correspond to which specific phonon mode..

Best Regards

shankar1729 commented 5 months ago

Please see the script within http://jdftx.org/EphMatrixElements.html that computes the Eliashberg spectral function on teh route to compute resistivity. In particular see this section of code from the penultimate block:

# Term combining matrix elements and velocity factors, summed over bands:
    w1 = np.exp(-0.5*((E1-mu)/Esigma)**2) * (1./(Esigma*np.sqrt(2*np.pi)))
    w2 = np.exp(-0.5*((E2-mu)/Esigma)**2) * (1./(Esigma*np.sqrt(2*np.pi)))
    term = np.einsum('ka,Kb,kKab,kKxab->kKx', w1, w2, vFactor, np.abs(g)**2)
    # Histogram by phonon frequency:
    F += np.histogram(omegaPh.flatten(), omegaPhBins, weights=term.flatten())[0]

In that term is calculated with indices kKx where x is phonon mode, and this is flattened and summed over in the next line. You can easily modify this section to separately histogram contributions from each phonon mode to get the mode-resolved spectral function.

Best, Shankar

dijasila commented 5 months ago

Hi Shankar, Apologies for my poor understanding of the code that calculates Eliashberg functions. Can you please elaborate a bit more on your response? So what you are saying is I should process the 'term' seperately for each phonon mode and not take a sum as is in line "F += np.histogram(omegaPh.flatten(), omegaPhBins, weights=term.flatten())[0]"?

something like this.... term = np.einsum('ka,Kb,kKab,kKxab->kKx', w1, w2, vFactor, np.abs(g)**2)

Histogram by phonon frequency:

F = np.histogram(omegaPh.flatten(), omegaPhBins, weights=term.flatten())[0] F[iblock] += F

shankar1729 commented 5 months ago

No, something like this (omitting the unmodified block of code in the middle, showing with ...):

F = np.zeros((nModes, len(omegaPhMid)))  # spectral function for each phonon mode
print('Sampling transport spectral function:', end=' ')
for iBlock in range(Nblocks):
    # Select blocks of k-points from above set:

    ...

    # Histogram by phonon frequency, separately for each phonon mode:
    for iMode in range(nModes):
        F[iMode] += np.histogram(omegaPh[..., iMode].flatten(), omegaPhBins, weights=term[..., iMode].flatten())[0]
    print(iBlock+1, end=' '); sys.stdout.flush()

where the additional dimension of F is each phonon mode. You have to of course modify the subsequent plotting and analysis routines accordingly. (I have not tested this explicitly; just sketching out what would need to be done.)