<jupyter_start><jupyter_text>a-Si 512 atom, getting eigenvector and frequency based on phonopy, then evaluate dynamic structure factor<jupyter_code>import numpy as np
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
from ase.io import read

atoms=read('optimized.vasp',format='vasp')


unitcell = PhonopyAtoms(symbols=atoms.symbols,
                    cell=atoms.cell,
                    positions=atoms.positions)
phonon = Phonopy(unitcell,
                [[1, 0, 0], [0, 1, 0], [0, 0, 1]],
                primitive_matrix=[[1,0,0],[0,1,0],[0,0,1]],log_level=1)
import pyAF.data_parse
FC_file='FORCE_CONSTANTS_2ND'
natom=len(atoms.positions)
#read as ndarray, no mass scaling here
fc_phonopy=pyAF.data_parse.read_fc_phonopy_noscale(FC_file,natom)
fc_phonopy.shape
phonon.set_force_constants(fc_phonopy)
phonon.symmetrize_force_constants(show_drift=True)
symetrized_fc=phonon.force_constants
#getting eigenvector and freq at k=0

mesh = [1, 1, 1]
phonon.run_mesh(mesh,with_eigenvectors=True)
mesh_dict = phonon.get_mesh_dict()
frequencies = mesh_dict['frequencies']
eigenvectors =mesh_dict['eigenvectors']
frequencies.shape, eigenvectors.shape<jupyter_output><empty_output><jupyter_text>phonopy yaml format: Eigenvector at each phonon mode. Each eigenvector e of dynamical matrix is shown as sets of three complex values of each atom along the Cartesian axes in the primitive cell. The real and imaginary values correspond to the left and right, respectively. A set of eigenvectors comprising all bands at a q-point forms a unitary matrix obtained as the result of numpy.linalg.eigh, i.e., LAPACK of routine _heevd. Therefore eigenvectors correspond to the column vectors of the unitary matrix.Eigenvectors corresponding to phonopy yaml output are obtained as follows.```if eigvecs is not None:    for eigvecs_at_q in eigvecs:        for vec in eigvecs_at_q.T:            print vec```<jupyter_code>k=np.array([0.1,0.0,0.0])
dsf=np.zeros(natom*3)
for i,vec in enumerate(eigenvectors.T):
    for j in range(natom):
        nstart=3*j
        nend=3*(j+1)
        disp=np.squeeze(vec)[nstart:nend]
        pos=atoms.positions[j]
        phase=np.exp(1.0j*np.dot(k,pos))
        polarlization=np.dot(k,disp)
        factor=phase*polarlization
        dsf[i]=dsf[i]+factor*np.conjugate(factor)


dsf[0:30]
nmodes=natom*3
dyn_file='Dyn.form'
lammps_dyn=np.loadtxt(dyn_file).reshape((nmodes,nmodes))
eigenvalue, eigenvector=np.linalg.eigh(lammps_dyn)
k=np.array([0.01,0.0,0.0])
dsf=np.zeros(natom*3)
for i,vec in enumerate(eigenvector.T):
    for j in range(natom):
        nstart=3*j
        nend=3*(j+1)
        disp=np.squeeze(vec)[nstart:nend]
        pos=atoms.positions[j]
        phase=np.exp(1.0j*np.dot(k,pos))
        polarlization=np.dot(k,disp)
        factor=phase*polarlization
        dsf[i]=dsf[i]+factor*np.conjugate(factor)
dsf[0:100]
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
bins=np.linspace(0,100,101)
bins
def broadening(omega1, omega2, sigma):
    delta=1.0/np.pi*sigma/((omega1-omega2)**2+sigma*sigma)
    return delta
sigma=1.0
smearing=np.zeros(len(bins))
for i,e in enumerate(bins):
    val=0.0
    for j,w in enumerate(omega):
        val=val+broadening(w,e,sigma)*dsf[j]
    smearing[i]=val
smearing
import matplotlib.pyplot as plt 
plt.plot(bins,smearing)<jupyter_output><empty_output><jupyter_text>志賀さんのコードをそのまま移植してみる<jupyter_code>a=atoms.cell[0,0]
np.pi*2/a
nom=100
nwn=10

ommin=np.min(omega)
ommax=np.max(omega)
om=np.linspace(ommin,ommax,nom)
Gmin=0.02
Gmax=10*Gmin
wn=np.linspace(Gmin,Gmax,nwn)

dsf_L=np.zeros([nwn,nom])
dsf_T=np.zeros([nwn,nom])

EL=np.zeros([nwn,nmodes])

ei=eigenvector.reshape([3,natom,nmodes])

q=np.array([np.ones(natom),np.zeros(natom),np.zeros(natom)])

q.shape
for ind, w in enumerate(omega):
    for ik in range(nwn):
        k=wn[ik]*q
        phase=np.exp(1.0j*np.dot(k,atoms.positions))
        EL[ik,ind]=EL[ik,ind]+np.linalg.norm(np.dot(np.dot(k,ei[:,:,ind].T),phase))**2/3

sigma=1.0
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)
dsf_L[0,:]
for ik in range(nwn):
    plt.plot(dsf_L[ik,:],om,label=str(ik))
plt.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1374128919860627, 'fraction_numerical': 0.02221254355400697, 'mean_word_length': 4.574029126213592, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '141734', 'license': 'MIT', 'n_tokens_mistral': 1917, 'n_tokens_neox': 1811, 'n_words': 328, 'path': 'aSi512_test/postprocess.ipynb', 'repo_name': 'eminamitani/thermal_conductivity_code', 'size': '4592'}