PWSCF计算晶体的声子色散曲线和声子态密度及其后处理步骤
已有 5950 次阅读
|个人分类:|系统分类:
，产生)；得到得到。【为什么要设置成小的，多小算小？】计算出来的每一个q点的本征频率可按上一个blog中提到的方法处理一下后画图。This example illustrates how to use pw.x and ph.x to calculate phonon frequencies at Gamma and X for Si and C in the diamond structure and &for fcc-Ni.The calculation proceeds as follows (for the meaning of the cited input variables see the appropriate INPUT_* file)1) make a self-consistent calculation (input=si.scf.in, output=si.scf.out)2) make a phonon calculation for gamma point (input=si.phG.in,output=si.phG.out). Note that outdir and prefix have the same & values as in the input for scf calculation. & Note that the q-point (the Gamma point in this case) is read & after the namelist inputph.
& Dynamical matrix will be saved for further analysis in fildyn (whose &default name is 'matdyn', here overwritten to si.dynG) & The Gamma point is the ONLY one for which the phonon code can be run &directly after the scf run (provided save file was defined). & For all other points in the BZ a non-scf preparatory run is needed, as it is shown here for the X point. 【对于其他非gamma点还是要进行非自洽计算】3) make a phonon calculation for X point (input=si.phX.in, output=si.phX.out). Note the option &lnscf=.true.&: this instructs & the code to perform a needed non-scf calculation first. & Note that the q-point (the X point in this case) is read after the namelist inputph. & Dynamical matrices will be saved for further analysis in fildyn &(whose default name is 'matdyn', here overwritten to si.dynX).
& The fildyn file contains the dynamical matrices for q-points that are equivalent by symmetry to the q-point given in input (in this case the three X points), therefore when calculating phonons in the & whole BZ the code needs to be run only on inequivalent points. & The phonon code can perform also a 'single-mode' calculation where & only a single atom is moved in a single direction and the corresponding column of the dynamical matrix is calculated. This may be useful when studying vibrational properties of defects and/or surfaces.【如何理解】 & This feature is illustrated in the following part of the test.4) Perform again a scf calculation for Si. Why? because the phonon &calculation performed in the previous step has overwritten the &original wavefunctions. This step can be avoided by saving the &content of outdir and restoring it after step 3)5) make a phonon calculation at X point for a single mode (input=si.phXsingle.in, output=si.phXsingle.out).
& Apart from the specification of the mode to be calculated & (modenum=3), the input is identical to case 4). The only difference &is that fildyn is not given a name, so the default one (matdyn) is &used. This is done in order not to overwrite the previous files.
& The total dynamical matrix is not computed and matdyn is actually & left empty.# self-consistent calculation for Ni with US-PPcat & ni.scf.in && EOF &control & &calculation='scf' & &restart_mode='from_scratch', & &tprnfor = .true. & &prefix='ni', & &pseudo_dir = '$PSEUDO_DIR/', & &outdir='$TMP_DIR/' / &system & &ibrav=2, celldm(1) =6.65, nat= &1, ntyp= 1, & &nspin=2, & &starting_magnetization(1)=0.5, & &degauss=0.02, & &smearing='mp', & &occupations='smearing', & &ecutwfc =27.0 & &ecutrho =300.0 / &electrons & &conv_thr = &1.0d-8 & &mixing_beta = 0.7 /ATOMIC_SPECIES Ni &58.6934 &Ni.pbe-nd-rrkjus.UPFATOMIC_POSITIONS Ni 0.00 0.00 0.00K_POINTS AUTOMATIC4 4 4 1 1 1EOF$ECHO & &running the scf calculation for Ni...\c&$PW_COMMAND & ni.scf.in & ni.scf.outcheck_failure $?$ECHO & done&# phonon calculation at X &【这是非自洽计算？】cat & ni.phX.in && EOFphonons of Ni at X &inputph &tr2_ph=1.0d-14, &prefix='ni', &amass(1)=58.6934, &fildyn='nix.dyn', &outdir='$TMP_DIR/', /0.0 0.0 1.0EOF$ECHO & &running the phonon calculation at X for Ni...\c&$PH_COMMAND & ni.phX.in & ni.phX.outcheck_failure $?$ECHO & done&# clean TMP_DIR$ECHO & &cleaning $TMP_DIR...\c&rm -rf $TMP_DIR/*$ECHO & done&$ECHO$ECHO &$EXAMPLE_DIR: done&备注：（1）动力学矩阵文件*.dyn0的作用：存储需要计算的q点数（2）声子谱q点选择很重要，选择的不好，分析计算会出现q not allowed；这种情况下，可取成一样的。一般依据：Pdos计算前 asr选项去掉（4）声子谱计算前自洽计算的目的，还是为了得到电荷密度 & 新版pwscf不需要进行非自洽【非自洽计算的目的，是为了得到费米面附近更精确的电荷密度值】；只有在计算电声耦合常数时，才要求加入非自洽计算 （5）PW计算声子谱，测试gama点很重：即用不同kpoints来测gaman点….无论采用phonopy还是pw计算声子谱前一定要在优化的过程中加入力的收敛！！！单纯的结构预测没有必要设置力的收敛….（6）计算声子谱时q点如何选取，是否也应该进行测试？为什么要设置q点？&[20:00:20] xiaoqiugood: 用446的k点做一下scf然后马上杀掉，out文件里面就有产生的k点的个数了（个数不是简单的等于4x4x6吗？）[20:01:25] lv Jian: 不等于，考虑对称性以后有很多K点事等价的，所以实际计算的要比4*4*6小[20:03:37] xiaoqiugood: 哦有道理，我说的是声子谱中的q点设置，你的意思是说：可以将设置好的q点当成k点，放入自洽计算文件中，跑一下，得出其实际数目，对吗？[20:04:34] lv Jian: 对，q点和k点事一样的，这么跑一下，就是想看看446的格子可以产生多少个点[20:05:24] xiaoqiugood: 有道理，我跑一下......[20:05:46] lv Jian: 然后调整格子，使实际计算的点大概20-30左右[20:07:30] xiaoqiugood: 好的，谢谢...这样的设置大致可以产生多少个.dyn文件？.dyn是一个原子产生一个吗？[20:09:34] lv Jian: 就是你产生的q点的个数[20:10:32] lv Jian: dyn文件的个数和你产生的q点的个数是一样的[20:10:59] xiaoqiugood: 哦，明白了.... (y)[20:11:22] lv Jian: (handshake)[20:11:53] xiaoqiugood: 34个k点可以接受哈[20:15:16] lv Jian: 稍多一点，如果计算快的还行，如果计算慢的话你就在改下网格Fcc 4 4 2产生8个k点 grep number of k points *out*number of k points= & & 8 gaussian broad. (Ry)= &0.0200 & & ngauss = &1 & & & & & & & & & & & cart. coord. in units 2pi/a_0 & & & k( & &1) = ( & 0..0000000 & 0.0000000), wk = & 0.0625000 & & & k( & &2) = ( & 0..5000000 & 0.5000000), wk = & 0.2500000 & & & k( & &3) = ( & 0..2500000 & 0.2500000), wk = & 0.2500000 & & & k( & &4) = ( & 0..2500000 & 0.7500000), wk = & 0.7500000 & & & k( & &5) = ( & 0..0000000 & 0.0000000), wk = & 0.1875000 & & & k( & &6) = ( & 0..0000000 & 0.5000000), wk = & 0.1250000 & & & k( & &7) = ( & 0..5000000 & 1.0000000), wk = & 0.2500000 & & & k( & &8) = ( & 0..0000000 & 0.5000000), wk = & 0.1250000&Q值设定网格，按该设定网格计算完声子谱后，采用tools包中的小工具处理，得到高对称性点和声子谱。a)klist.ini &文件该文件，可通过xcrysden选点产生，文件格式如下：b)./cry-to-car.x &klist.ini→klist.out转换成晶体坐标c)./kpath-30.x &在高对称性点中插入的点数？d)./kband.x在klist.out中找出高对称性点，并标注出来&e)将klist.out中点放入phon-spectru中….执行后得到所要处理的声子谱文件1）计算得到.dyn文件；pw.x自洽，ph.x声子，得到的dyn文件数量等于产生的q点数量2）以下都是后处理：3）准备高对称性k点，按格式写入klist.ini文件4）执行cry-to-car.x程序，将klist.ini文件中的k点坐标由晶体坐标转换为笛卡尔坐标，得到相应的klist.out文件5）在各高对称性k点间插入间隔k点数：即，执行 ./kpath_30.x 6）运行./kband.x，在各k点加加入逗号（供后续程序读取数据方便）7）将klist.out文件内容放入phonon_spectra2.1.5文件中，按格式修改相关参数，执行程序，得到…. &运行，得到可用于绘图的文件【注意采用该程序提取的频率的单位为，不是】1. pw.x 和ph.x脚本自己先选定一些高对称点，并产生这些高对称点之间其他点。2. & 后处理步骤脚本&cat & matdyn.in&& EOF &input & &asr='simple', &amass(1)=1.008, amass(2)=44.955910 & &flfrc='a225.fc', flfrq='a225.freq' /151 &-1.000000 & & &0.5000000 & & &0. & & & & & 1 -0.9833333 & & &0.5000000 & & &1. & & & & & 2 -0.9666667 & & &0.5000000 & & &3. & & & & & 3 .......(略) -0.7666667 & & &0.7333333 & & &0. & & & & 149 -0.7583333 & & &0.7416667 & & &0. & & & & 150 -0.7500000 & & &0.7500000 & & &0. & & & & 151EOF $PW_ROOT/matdyn.x & matdyn.in& matdyn.out#这么多k点来自产生k点的小脚本# flfrc：用来设置输出力常数矩阵的文件; & &# flfrq：用来输出频率的文件#这里要输入151个特殊q点的坐标。与计算能带结构时一样，需先选出要计算的高对称q点的走向以及高对称点的坐标，然后产生这些线上的q点的坐标。计算出来的每一个q点的本征频率可按上一个blog中提到的方法处理一下后画图。#声子谱q点选择很重要，选择的不好，分析计算会出现q not allowed；这种情况下，可取成一样的。一般依据：pdos计算前 asr选项去掉；声子谱计算前自洽计算的目的，还是为了得到电荷密度。新版pwscf不需要进行非自洽【非自洽计算的目的，是为了得到费米面附近更精确的电荷密度值】；只有在计算电声耦合常数时，才要求加入非自洽计算。# dynmat.x applies various kinds of Acoustic SumRule (ASR), calculates LO-TO splitting at q = 0 in insulators, IR and Raman crosssections(if the coefficients have been properlycalculated), from the dynamical matrix produced by ph.x. &Notice again that asr='simple' makes the frequency ofacoustic phonons vanish but has no effect on the other modes&&&cat & phdos.in&&EOF &input & &asr='simple', & &dos=.true., & &amass(1)=1.008, & &amass(2)=44.955910, & &flfrc='a225.fc', & &nk1=20,nk2=10,nk3=10, & &deltaE=1.0, & &fldos='a225.dos'ndos=50 /EOF $PW_ROOT/matdyn.x & phdos.in& phdos.out&#要计算声子态密度,dos必须设置为.true.另外fldos用来设置输出的态密度值(文件)，计算态密度时要用更密的q点网格，这需设置nk1, nk2, nk3。另外还有态密度的能量刻度上的点的数目，可由ndos来设置(非必须)。# fldos：用来存储声子态密度的文件;表示例子文件：# job script for quantum espresso# the name of queue, should be cscore or csnodeAPP_NAME=&cscore&# number of processesNP=16# fix command styleRUN=&RAW&# current directoryCURDIR=$PWD# build .nodelistrm -rf $CURDIR/.nodelist && /dev/nullfor i in `echo $LSB_HOSTS`doecho $i && $CURDIR/.nodelistdone# directory for quantum espressoEXEDIR=$HOME/bin/espresso502/bin# directory for mpi (mvapich)MPIDIR=/home/compiler/mpi/mvapich/1.0/icc.ifort-11.1/bin# start calculation (scf)#$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/pw.x -in scf.in & scf.out# start calculation (nscf)#$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/pw.x -in nscf.in & nscf.out# start calculation (phonon)#$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/pw.x -in ph.in & ph.out#PARA_PREFIX='$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist'#export PARA_PREFIX####################################################################PW_ROOT=/home/users/spclibtt/bin/espresso502/binPSEUDO_DIR=/home/users/spclibtt/home/yexq/shanghai-opt/pw-work/pseudoTMP_DIR=/home/users/spclibtt/home/yexq/shanghai-opt/pw-work/environment/tmp#export PARA_PREFIX='mpirun -np 2'# or & export PARA_PREFIX='mpirun' ,export &PARA_POSTFIX= -np 3export PATH=/home/compiler/mpi/mvapich/1.0/icc.ifort-11.1/bin:$PATHexport PARA_PREFIX='mpirun -np 16'export PW_ROOT PSEUDO_DIR TMP_DIRexport name='sch4-i4mm'rm -r $TMP_DIR/*.save rm $TMP_DIR/*###########self-consistent calculation ######################for a in 2500docat & $name.scf.in_$a && EOF &control & &calculation = 'scf' & &restart_mode='from_scratch', & &prefix='$name', & &pseudo_dir = '$PSEUDO_DIR/', & &outdir='$TMP_DIR/' & &tstress=.t., & &tprnfor=.t./&system &ibrav=0, &nat=10, &ntyp=2, &ecutwfc=90.0, ecutrho = 760.0 &nbnd= 30, &occupations ='smearing', degauss =0.01 &smearing ='mp'/ &electrons & &mixing_beta = 0.7 & &conv_thr = &1.0d-8 /ATOMIC_SPECIESH & & 1.008 & & &H.pbe-van_ak.UPFSc & &44.9559 & &Sc.pbe-nsp-van.UPFCELL_PARAMETERS {bohr}4. & & -4.27115E-08 & &0-4.27115E-08 & &4. & & -9.49146E-094.74573E-09 & & -1.89829E-08 & &8.ATOMIC_POSITIONS {crystal}H & & & &0. & 0. & 0.H & & & &0. & 0. & 0.H & & & -0. & 0. & 0.H & & & -0. & 0. & 0.H & & & -0. & 0. & 0.H & & & -0. & 0. & 0.H & & & &0. & 0. & 0.H & & & &0. & 0. & 0.Sc & & &-0. &-0. & 0.Sc & & & 0. & 0. & 0.K_POINTS {automatic}16 16 9 0 0 0EOF$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/pw.x && $name.scf.in_$a &$name.scf.out_$adonecat &$name.ph.in && EOFphonon of $name &at &X&inputphtr2_ph=1.0d-12,prefix='$name',fildvscf='$name.dv',amass(1)=1.008amass(2)=44.955910,trans=.true.,ldisp=.true.,outdir='$TMP_DIR',fildyn='$name.dyn',nq1=4,nq2=4,nq3=2 /EOF$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/ph.x && $name.ph.in & $name.ph.out完整的后处理程序#!/bin/bash####################################################################. ../../../../../environment/environment_1# export PARA_PREFIX='mpirun &-np 6 '# or & export PARA_PREFIX='mpirun' ,export &PARA_POSTFIX= -np 3export &TMP_DIR=../../../../../tmp/tmp5 rm -r $TMP_DIR/*.save rm $TMP_DIR/*export name='c225'#############色散曲线#############cat & q2r.in && EOF &input & fildyn='c225.dyn', zasr='simple', flfrc='c225.fc' /EOF $PW_ROOT/q2r.x & q2r.in & q2r.outcat & matdyn.in && EOF &input & &asr='simple', &amass(1)=1.008,amass(2)=44.955910 & &flfrc='c225.fc', flfrq='c225.freq' /151 &-1.000000 & & &0.5000000 & & &0. & & & & & 1 -0.9833333 & & &0.5000000 & & &1. & & & & & 2 -0.9666667 & & &0.5000000 & & &3. & & & & & 3 ....... -0.7500000 & & &0.7500000 & & &0. & & & & 151EOF $PW_ROOT/matdyn.x & matdyn.in & matdyn.out###################dos##########################cat & phdos.in &&EOF &input & &asr=.true., & &dos=.true., & &amass(1)=1.008 & &amass(2)=44.955910, & &flfrc='c225.fc', & &nk1=20,nk2=10,nk3=10, & &deltaE=1.0, & &fldos='c225.dos' /EOF $PW_ROOT/matdyn.x & phdos.in & phdos.out
转载本文请联系原作者获取授权，同时请注明本文来自叶小球科学网博客。链接地址：
上一篇：下一篇：
当前推荐数：1
评论 ( 个评论)
扫一扫，分享此博文
作者的精选博文
作者的其他最新博文
热门博文导读
Powered by
Copyright &小木虫 --- 500万硕博科研人员喜爱的学术科研平台
&&查看话题
Quantum-ESPRESSO 5.0.3，Windows 32-bit，串行版
:D 在MinGW下，成功使用gfortran和gcc编译了串行版的PW 5.0.3（已经打上了5.0.2 to 5.0.3的补丁）。
当然，由于又是Windows，又是串行，又是GNU编译器，不要对效率有太高的要求了。事实上，指望用这个来做计算工作基本不太可能。不过，用这个版本在Linux编译之前学习一下PWscf的使用，试用一下Linux下软件使用的感受，还是不错的，也可以做到软件学习和编译学习两手抓嘛。
本次提供的编译后的模块分别有：
1，PWscf的基本程序：pw.exe
2，声子计算模块
3，能带计算模块
4，态密度计算模块
5，电荷密度绘制模块
6，NEB模块
7，TDDFPT模块
8，赝势生成模块
顺便附加了一个小例子，计算AlAs的能带结构，可以按照：
pw.exe & AlAs.scf.in
pw.exe & AlAs.nscf.in
bands.exe & AlAs.bands.in
plotband.exe & AlAs.plotband.in
的顺序进行计算。
所谓的“安装”方法与帖：http://emuch.net/bbs/viewthread.php?tid=4923069&fpage=3 中所述相同，请那些只喜欢双击却从不看文献/教程的懒人们不要硬生生来句“双击无法运行”，把自己的全责推成他人的错误。
顺便多说一句，既然是做计算，基本的计算机知识总要自己补吧？就好像做电镜的肯定要知道电镜基本结构和机理，开汽车的在驾校理论课上还得讲解一点汽车的基本知识，做计算却不去了解计算机，连自己的吃饭的工具（至少是研究生期间）都懒得关心，是不是太说不过去了呢？
计算后的AlAs能带结构图：
现在Quantum Espresso的官方下载网站提供了最新版QE的Windows版本，分别为32bit和64bit的串行及并行版。对于本贴中的资源可以无视了。
:D 从来没用过NEB，要不hakuna测试一下呗？
第一个问题可以看黄良峰的blog及他提供的小工具。
第二个问题，本论坛讨论很多，可以自行搜索一下。
谢谢您&&我看到他分享的一个将声子谱按对称性分配的脚本
谢谢&&我找到了
一直想用一下QE，但一直没有时间动手去弄这个。
我只是根据VASP使用NEB时必须并行的情况有上面一问，从没用过QE, souledge兄高抬了
第一个问题可以看黄良峰的blog及他提供的小工具。
能提供一下网址吗？谢谢！
并行编译时无法在MinGW下跑的，何况，最新的QE 5.1.2已经官方提供了32和64位的Windows版，有串行也有并行版。
研究生必备与500万研究生在线互动！
扫描下载送金币
浏览器进程
打开微信扫一扫
随时随地聊科研& 记住密码
1) k点和q点的物理意义；
2) k点和q点在倒空间的取样方法；
3) 电声耦合作用计算时，k点和q点的取点要求
1.nq,xq的区别
nq属于放入这个 INPUTPH 字段中的参数
xq(1) xq(2) xq(3)
The phonon wavevector, in units of 2pi/a0
(a0 = lattic ...
1） 零点振动能的本质；
Perosonal view:
ZPE是物质的固有属性，不随计算方法而改变，即选取不同的参数如赝势种类、截断能大小等，只要精度足够，不影响计算得到的零点能数值。
2） 零点能的表达式；
3） 从PWscf获取零点能的步骤
感谢LHY分享此脚本：
#!/usr/bin/python
7.QE中的对称性运用与自洽、声子计算
Author: paulatto
14:32:05 +0100 (Tue, 20 Nov 2012)
New Revision: 9630
trunk/espresso/PHonon/PH/set_irr_sym.f90
trunk/espresso/PHonon/PH/sgam_ph.f90
Short explanation: minus_q symmetry fixed, for many years
1） 线性响应计算声子谱的原理；
2） gamma点测试的重要性及gamma点与声学支的关系；
3） PWscf线性响应计算声子谱的步骤
由此可见gamma点测试的结果是快速判定结构是否有虚频，如有可与光学支最低频率数值可比拟的虚频
则可认为结构不稳定！！！
自洽k点取20 20 11
版计算电声耦合常数分为以下几个步骤：
1. make a self-consistent calculation for Alusing
a dense grid of k-points.
The dense grid must contain all kand k+q grid points used in the subsequent electron-phonon calculation and must be denseenough to
produce accurate
转载：PWscf
谱的步骤和注意事项
、建立你所需要的晶体结构
在这步中，最重要的是测试
是否能找到晶体的全部对称性；
等，这取决于你是否想让你计算的结构加压力，或者应力什么的。
前面两步是前提，但 ...
1） 自洽计算的目的是为了check结构建立是否合理，程序能否顺利运行；
2） kpoints测试的重要性； 选择合适的参数，可在保证质量的情况下，节省计算时间
3) ecut 测试的重要行
经验总结：
1） 赝势选取后，测试赝势的截断能，最初值设置参考依据是…..？从赝势文件中能看 ...
优化计算例子文件
--run_optim
####################################################################
./home/xiaoqiugood/pwscf/pw/environment/environment_1
exportPARA_PREFIX='mpirun -np 6'
export PARA_PREFIX='mpirun' ...
请教一下关于ecut、kpoint和degauss的设定 （1）ecut和kpoint的设定需要收敛测试，是不是只要ecut和kpoint的选定使总能收敛了就行了呢？需不需要测试力的收敛呢？ （2）对于pwscf里面的ecut，ecutwfc 和ecutrho 需要同时测试呢，还是仅仅测试ecutwfc然后把ecutrho选定为ecutwfc的4倍以上呢？ （3）K_POINTS (automatic) & ...
1) 声子谱的计算原理；
2) PWSCF计算声子谱的步骤；
3) 高对称性点的插入及后处理步骤。
声子谱计算步骤小结：以下两个步骤写入
脚本文件一次执行
进行电子密度的自洽计算
网格点进行动力学矩阵元的计算
fildyn=&# ...
Powered by
Copyright &[转载]PWscf计算IR和Raman谱的步骤和注意事项
1、建立你所需要的晶体结构
在这步中，最重要的是测试
2）pw.x是否能找到晶体的全部对称性；
2、vc-relax&
&relax等，这取决于你是否想让你计算的结构加压力，或者应力什么的。
前面两步是前提，但这里就不多说了，我主要谈谈如何计算红外和拉曼谱
3、用优化好的晶体结构进行scf计算
&&&title =
&& & prefix =
&&&calculation =
restart_mode = 'from_scratch' ,
&& &wf_collect =
&& & outdir =
&& &pseudo_dir =
&etot_conv_thr = 1.0d-5 ,
&forc_conv_thr = 1.0d-4 ,
&& &tstress =
&& &tprnfor =
&&&nstep = 100
&&&ibrav =
&& & celldm(1) =
&& & celldm(3) =
&& & nat =
&& &ntyp =
&& &ecutwfc =
&&&nosym =
&ELECTRONS
&electron_maxstep = 100,
&&&conv_thr =
ATOMIC_SPECIES
&16-S.GGA.fhi.UPF
Zn&&65.409&
&30-Zn.GGA.fhi.UPF
ATOMIC_POSITIONS crystal
&-0.&&-0.&
K_POINTS automatic
&&16 16 16 0 0 0
结果文件中，有
& & 12 Sym.Ops. (no
inversion)
这说明pw.x找到了12个对称操作。如何判断pw.x找到了全部对称操作呢？我是根据materials
studio判断的，这个问题有空再讨论。
总之，pw.x一定要找到全部的对称操作才行。
除了要找到全部对称性外，另外就是使用较密的k-mesh，为后续的声子性质计算做准备；
还有一个需要注意的，我这里用了模守恒赝势，是从abinit网站上下的，然后用fhi2upf.x工具转换的，在pwtools下有。
4、用ph.x计算gamma点的声子性质，即一阶红外和拉曼性质
ph.x & ph.in & ph.out
ph.in 文件如下
phonon calculation for ZnS
&&tr2_ph=1.0d-12,
&&prefix='ph',
&&epsil=.true.,
&&trans=.true.,
&&lraman=.true.,
&&amass(1)=32.066,
&&amass(2)=65.409,
&&outdir='./',
&&fildyn='ZnS_IR.dynG',
&&fildrho='ZnS_IR.drho',
0.0 0.0 0.0
这里我要说明下
&&tr2_ph=1.0d-12,
这好像是默认值，可以调高些，计算量增加的不多，但是能够提高声子的计算精确度；
&&epsil=.true.,
&&trans=.true.,
&&lraman=.true.,
&&这三个设置成true，就能计算红外和拉曼了
&&fildyn='ZnS_IR.dynG',
这个是输出的动力学矩阵文件，在后处理中需要对该文件处理，大家可以看看里面是什么内容；
下面贴出ph.out 文件的最后输出信息
**************************************************************************
&&&omega( 1)
& 1.007969 [THz] =&
&& &33.622437
&&&omega( 2)
& 1.274035 [THz] =&
&& &42.497508
&&&omega( 3)
& 1.274035 [THz] =&
&& &42.497508
&&&omega( 4)
& 2.280533 [THz] =&
&& &76.070895
&&&omega( 5)
& 2.280533 [THz] =&
&& &76.070895
&&&omega( 6)
& 6.220740 [THz] =&
&&&207.502947
&&&omega( 7)
& 8.545154 [THz] =&
&&&285.037559
&&&omega( 8)
& 8.616679 [THz] =&
&&&287.423368
&&&omega( 9)
& 8.616679 [THz] =&
&&&287.423368
&&&omega(10)
& 8.870194 [THz] =&
&&&295.879779
&&&omega(11)
& 8.870194 [THz] =&
&&&295.879779
&&&omega(12)
&10.242682 [THz] =&
&&&341.661368
**************************************************************************
symmetry, C_6v (6mm)&&point
&&&omega(&&1
&33.6&&[cm-1]&
&&&omega(&&2
&42.5&&[cm-1]&
&&&omega(&&4
&76.1&&[cm-1]&
&&&omega(&&6
&&&207.5&&[cm-1]&
&&&omega(&&7
&&&285.0&&[cm-1]&
&&&omega(&&8
&&&287.4&&[cm-1]&
&&&omega( 10 -
&&&295.9&&[cm-1]&
&&&omega( 12 -
&&&341.7&&[cm-1]&
**************************************************************************
频率大小，对称性都给出了，不可约表示也给出来了
但是红外，拉曼的强度没有给出
下面进行下一步！
5，用dynmat.x 计算红外和拉曼强度
dynmat.x & ZnS001.in &
ZnS001.out
ZnS001.in文件
fildyn='ZnS_IR.dynG',
asr='crystal',
q(1)=0,q(2)=0,q(3)=1,
filout='ZnS_001.out'
filmol='ZnS_001.mold'
filxsf='ZnS_001.axsf'/
计算得到的
ZnS001.out文件
&&&IR cross
sections are in (D/A)^2/amu units
&&&Raman cross
sections are in A^4/amu units
&&&multiply
Raman by 0.140845 for Clausius-Mossotti correction
&-0.0000& &
0.0000& &&
&0.0155& &
0.0000& &&
&0.0372& &
0.0000& &&
&0.0132& &
& 2.2805& &
0.0000& &&
&0.1087& &
& 2.2805& &
0.0000& &&
&0.1087& &
& 207.50& &
0.0000& &&
&0.0000& &
& 288.35& &
8.1445& &&
&1.6782& &
& 288.35& &
8.1445& &&
&1.6782& &
& 295.88& &
0.0000& &&
&2.3796& &
& 295.88& &
0.0000& &&
&2.3796& &
&10.2427& &
0.0000& &&
&0.0000& &
&10.5003& &
9.0605& &&
&3.2269& &
IR 和 Raman 栏里就是表示计算的强度
可以自己写个小的进行高斯展宽的程序，那样就能画个拉曼谱了
如果有需要这个小程序的，可以站内消息联系。这里就不多谈这个小程序了
我的感觉：
1，计算红外和拉曼对系统的要求比较高，不是什么系统都是适合用来计算红外拉曼性质的。
比如强关联的体系，因为不好找赝势，而且需要很多测试。
2，红外和拉曼的计算量还是挺大的，ph.x的运行时间比较久，所以大的系统，比如原子数比较多，或者对称性很低的晶体结构，
可能需要的时间更多。
3，小分子的红外和拉曼我没有计算过，但我觉得还是可以计算的。不过有些特殊的结构，比如有5重对称性的分子，C60，估计不好
计算。因为晶体群对称性中没有5重对称性，这样的话，由于对称性的减小，计算红外和拉曼将要花太多的时间。如果有计算过小分子
红外和拉曼的虫友，可以互相交流下。
4，提供下我发表的几篇简单的工作和我经常看的和红外拉曼计算相关的文献
& Y. C. Cheng, H. T. Chen, X. X. Li, X. L. Wu, J.
Zhu, S. H. Li, and P. K. Chu, J. Appl. Phys. 105, 083511
& Y. C. Cheng, C. Q. Jin, F. Gao, X. L. Wu, W.
Zhong, S. H. Li, and P. K. Chu, J. Appl. Phys. 106, 123505
& J. Sun, X. F. Zhou, J. Chen, Y. X. Fan, H. T.
Wang, X. J. Guo, J. L. He, and Y. J. Tian, Phys. Rev. B 74, 193101
& M. Lazzeri, and F. Mauri, Phys. Rev. Lett. 90,
& S. Baroni, S. de Gironcoli, A. Dal Corso, and P.
Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
以上网友发言只代表其个人观点，不代表新浪网的观点或立场。