<jupyter_start><jupyter_text>Dynamical Structure Factor(Still someting wrong)<jupyter_code>from ase.io import read

atoms=read('optimized.vasp',format='vasp')
import numpy as np
natom=len(atoms.positions)
nmodes=natom*3
dyn_file='Dyn.form'
lammps_dyn=np.loadtxt(dyn_file).reshape((nmodes,nmodes))
eigenvalue, eigenvector=np.linalg.eigh(lammps_dyn)
def broadening(omega1, omega2, sigma):
    delta=1.0/np.pi*sigma/((omega1-omega2)**2+sigma*sigma)
    return delta
a_super=atoms.cell[0,0]
a0=a_super/4
a0
kmin=2.0*np.pi/21.8893169999999984
kmin
omega_min=5.0
omega_max=100.0
nom=100
nwn=10
om=np.linspace(omega_min,omega_max,nom)
wn=np.linspace(kmin,kmin*10,nwn)

dsf_L=np.zeros([nwn,nom])
EL=np.zeros([nwn,nmodes])
wn
import pyAF.constants
pc=pyAF.constants.physical_constants()
omega=[]
#extract minimum index of negative frequency
mode_negative=0
for i in range(nmodes):
    if eigenvalue[i] <0.0:
        val=-np.sqrt(-eigenvalue[i])*pc.scale_cm
        omega.append(val)
        mode_negative=i
    else:
        val=np.sqrt(eigenvalue[i])*pc.scale_cm
        omega.append(val)
omega[100]
k=np.array([[1,0,0]])
np.repeat(k,3,axis=0)
#test for single mode eigenvector
ind=10
partial=eigenvector[:,ind]
print(partial)
EL=np.zeros([nwn,nmodes])

for ik in range(nwn):
    #100 direction
    k=np.array([[0,wn[ik],0]])
    n_1=np.array([[1.0,0.0,0.0]])
    disp=np.reshape(partial,(natom,3))
    polarlization=np.dot(n_1,disp.T)
    phase=np.exp(1.0j*np.dot(k,atoms.positions.T))
    value=np.dot(polarlization,phase.T)
    
    EL[ik,ind]=EL[ik,ind]+np.abs(value)**2/3
EL[:,ind]
import matplotlib.pyplot as plt 
plt.plot(wn,EL[:,ind])
omega[ind]
b=[]
sigma=0.01
for iw in range(nom):
    b.append(broadening(om[iw],omega[ind],sigma))
b
dsf_L=np.zeros([nwn,nom])
sigma=0.5
for iw in range(nom):
    for ik in range(nwn):
        dsf_L[ik,iw]=dsf_L[ik,iw]+EL[ik,ind]*broadening(om[iw],omega[ind],sigma)
dsf_L[0,:]
import matplotlib.pyplot as plt 
for ik in range(nwn):
    plt.plot(dsf_L[ik,:],om,label=str(ik))
plt.legend()
#[100] direction
dsf_L=np.zeros([nwn,nom])
EL=np.zeros([nwn,nmodes])

for ind, vec in enumerate(eigenvector.T):
    for ik in range(nwn):
        k=np.array([wn[ik],0,0])
        contrib=0.0
        for j in range(natom):
            nstart=3*j
            nend=3*(j+1)
            disp=vec[nstart:nend]
            pos=atoms.positions[j]
            phase=np.exp(1.0j*np.dot(k,pos))
            polarlization=np.dot(k,disp)
            factor=phase*polarlization
            contrib=contrib+factor

        EL[ik,ind]=EL[ik,ind]+np.abs(contrib)**2/3

EL[0,:]
sigma=0.01
for iw in range(nom):
    for ik in range(nwn):
        for ind in range(nmodes):
            dsf_L[ik,iw]=dsf_L[ik,iw]+EL[ik,ind]*broadening(om[iw],omega[ind],sigma)
np.linalg.norm(k)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1737891737891738, 'fraction_numerical': 0.03596866096866097, 'mean_word_length': 4.320075757575758, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 1, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '320166', 'license': 'MIT', 'n_tokens_mistral': 1259, 'n_tokens_neox': 1212, 'n_words': 170, 'path': 'aSi512_test/dsf.ipynb', 'repo_name': 'eminamitani/thermal_conductivity_code', 'size': '2808'}