ZZ:处理PWSCF中计算得到本征频率

转载自：http://valenhou.blog.edu.cn/2007/133199.html

处理PWSCF中计算得到本征频率　(2007-02-16 17:32)

• 　　分类：PWSCF　标签： pwscf   

在采用pwscf中的ph.x以及matdyn.x计算得到了沿高对称q点的本征频率后，需将数据转换一下，以画出声子色散曲线。下面是计算出来的本征频率文件中的格式：
&plot nbnd= 6, nks= 91 /
0.500000 0.000000 0.500000
139.6336 194.6691 277.0816 409.1945 460.9638 475.5246
0.475000 0.000000 0.475000
136.1839 189.6126 269.1646 418.0272 463.4473 475.7502
0.450000 0.000000 0.450000
132.4453 184.0388 260.9538 426.7985 466.0779 475.9910
0.425000 0.000000 0.425000
128.4070 177.8615 252.3483 435.3956 468.8332 476.3233
0.400000 0.000000 0.400000
124.0518 171.0092 243.2433 443.7172 471.6875 476.8231
0.375000 0.000000 0.375000
119.3569 163.4369 233.5371 451.6723 474.6123 477.5588
0.350000 0.000000 0.350000
114.2958 155.1312 223.1375 459.1796 477.5767 478.5845
0.325000 0.000000 0.325000
.......................
第一行中的nbnds告诉了有每个q点有多少个本征振动频率，nks告诉了计算了多少个q点。下面的接着的每两行是：
q点的坐标
每个q点对应的本征振动频率

要将q点的坐标转换为2D图中的横坐标上的长度值，下面的fortran90代码就是将这样的计算数据转换为q的长度值以及对应的本征频率值，以便画图。
! for PWSCF to plot phonon dispersion
!
program prog
real, allocatable :: e(:,:)
real, allocatable :: k(:,:)
real, dimension(3) ::k0,a
character(len=32):: input, output
character(len=32):: xx, yy

!read(5,*) nbands
!write(6,*) 'number of bands to be plotted'
!read(5,*) nk
write(6,*) 'name of input file, name of output file'
read(5,*) input,output

open(10,file=input, status='old')
open(11,file=output, status='new')

read(10,*) xx, yy, nbands, yy, nk

allocate(e(nk,nbands))
allocate(k(nk,3))
do i=1,nk
read(10,*)(k(i,j),j=1,3)
write(6,*)(k(i,j),j=1,3)
read(10,*) (e(i,n),n=1,nbands)
write(13,9030) (e(i,n),n=1,nbands)
9030 format (8f9.5)
enddo
do j=1,nbands
dk=0
do i=1,nk
if (i.eq.1) then
k0=k(i,:)
endif
a=k(i,:)-k0
dk=dk+sqrt(dot_product(a,a))

write(11,*)dk,e(i,j)
k0=k(i,:)
enddo
write(11,*)
enddo
stop
end program prog
--------
这是根据pwscf提供的一个代码所改的。