Support s(q,w) format and not just s(alpha,beta)

[tkittel]

Discussing with @marquezj , @dddijulio , @xushuqi77 . It seems S(q,w) = (1/kT)S(alpha,beta), and beta has opposite sign of omega (beta=-hbaromega/kT, alpha = hbar^2/2MkT * qval^2).

We need to specify units for q and omega.

[marquezj]

@tkittel we were discussing with @dddijulio that in some applications $S(q, \omega)$ can include the scattering lengths. This is something important to consider in conversion because usually $S(\alpha, \beta)$ is used in the incoherent approximation, and the scattering lengths or bound scattering cross sections are factored out.

[tkittel]

Good point, which also touches upon the single-phonon kernels promised by @xxcai1 . In that case, it might be better to keep the current support for s(alpha,beta) kernels in C++ and NCMAT for now. We could, however, provide some utilities for conversions to/from other formats in a python module. If anyone can write up some decent proposals for what functions might be useful, I'll be happy to include them in a future release.

[xxcai1]

I happen to have a function for that. The function is very very short, but the math is a bit long

I can imaging direct sampling sqw is basically writing a new class and can take a very long time to debug.

# Q^2 =  k_i^2 + k_f^2 -2*k_i*k_f*cos(\theta)
# so QdQ = k_f*k_i sin(\theta) d\theta --> sin(\theta) d\theta = QdQ/(k_f*k_i)

# given that d\Omega = sin(\theta) d\theta d \phi
# d\Omega = Q/(k_i*k_f) dQ d\phi

# totxs = \int \int k_f*Sqw/k_i d\Omega dE'
# -->  totxs = 2*pi* \int \int Q*Sqw/k_i/k_i dQ dE' , for isotropic scatters

# sqw2sab converts the Sqw to the Sab in the eq.4 in the rejection paper using the Jacobian, sab=sqw det(J) 
# NOTICE, THE POWER OF SCATTERING LENGTH IS BEING ABSORBED IN THE SQW

# Q in Aa^-1
# angular frequency in Hz

import numpy as pi
const_eV2kk = 1/2.072124652399821e-3
const_hbar = 6.582119569509068e-16
        
def sqw2sab(sqw, Q, fre, kt):
    alpha = Q*Q/(kt*const_eV2kk)
    beta = fre*const_hbar/kt
    sab = sqw*0.5*kt*kt*const_eV2kk/(2*np.pi*const_hbar)
    return sab, alpha, beta

[XuShuqi7]

Hi @xxcai1, I would like to ask a question concerning the conversion between sab and sqw in the sqw2sab function. If I well understood, we have: $\alpha=\dfrac{Q^2}{Ck_BT}$ and $\beta=\dfrac{E^\prime-E}{k_BT}=\dfrac{\hbar\omega}{k_BT}$ where $C$ is a conversion constant from eV to $\mathring{\mathrm{A}}^{-2}$. Then we can get: $\textrm{d}\alpha=\dfrac{2Q\textrm{d}Q}{Ck_BT}$ and $\textrm{d}\beta=\dfrac{\textrm{d}(\hbar\omega)}{k_BT}$. Therefore, the Jacobian determinant is: $\dfrac{C(k_BT)^2}{2}$. However, in the sqw2sab function, you multiplied by an extra term $\dfrac{1}{2\pi\hbar}$. I would like to know why we multiply this term, thank you!

[xxcai1]

I think that depends on the unit system that you used to generate the sqw. Likely, the planck number is unity in the code I'm using.

[XuShuqi7]

@xxcai1 Thank you for your answer! Indeed, if the Planck constant is set to be unity, it is consistent.

I would also like to ask what unit is conventionally used for $S(Q,\omega)$? If we follow the definition in the NCrystal paper, i.e, $\dfrac{\textrm{d}^2\sigma}{\textrm{d}\Omega\textrm{d}E^\prime}=\dfrac{k_f}{k_i}S(Q,\omega)$, $S(Q,\omega)$ is it in barn/eV/rad? I check in the OCLIMAX paper, but it seems that the arbitrary unit is used.

[xxcai1]

I don't know about the software. But in phonopy, you can see the units are different for different DFT codes. https://phonopy.github.io/phonopy/interfaces.html?highlight=unit. You may want to check if OCLIMAX is using the unit system of the DFT code or its own one.

[XuShuqi7]

Hi @xxcai1, thank you for your answer! It seems that the link cannot be opened, but I check in the documentation of Phonopy which mentions that the unit of $S(Q,\omega)$ is $\textrm{m}^2$/J (http://phonopy.github.io/phonopy/dynamic-structure-factor.html#dynamic-structure-factor).

Excuse me for my last reply which is a little ambiguous, in fact, OCLIMAX was mentioned just as an example. What I tried to do is to convert the $S(\alpha,\beta)$ generated by NCrystal to $S(Q,\omega)$. Since we have: $\dfrac{\textrm{d}^2\sigma}{\textrm{d}\Omega\textrm{d}E^\prime}=\sqrt{\dfrac{E^\prime}{E}}S(Q,\omega)$ and $\dfrac{\textrm{d}^2\sigma}{\textrm{d}\Omega\textrm{d}E^\prime}=\dfrac{\sigma_b}{4\pi k_BT}\sqrt{\dfrac{E^\prime}{E}}S(\alpha,\beta)$, the double differential cross section is the same no matter $S(Q,\omega)$ or $S(\alpha,\beta)$ is used, so we have the conversion: $S(Q,\omega)=\dfrac{\sigma_b}{4\pi k_BT}S(\alpha,\beta)$.

This was the conversion I used. However, the conversion mentioned by you uses the Jacobian determinant. The results of these two approaches are so different. I would like to know if you any idea about this? Thx!

[ramic-k]

@xxcai1 @XuShuqi7 Hi guys, if you read the OCLIMAX manual (https://docs.google.com/viewer?a=v&pid=sites&srcid=ZGVmYXVsdGRvbWFpbnxvcm5saWNlbWFufGd4OjZmMjBlY2IxYTQxYWRhMDQ) when NORM=1 flag is used (which is default) S(Q,W) is normalized to barn/energy (whatever your choice for energy units might be in OCLIMAX parameters file). Hopefully this helps a bit

[XuShuqi7]

Hi @ramic-k, thank you for the doc! I think I get the answer for the absolute unit of $S(Q,\omega)$.