1) 声子谱的计算原理；

2) PWSCF计算声子谱的步骤；

3) 高对称性点的插入及后处理步骤。

声子谱计算步骤小结：以下两个步骤写入prun_phon脚本文件一次执行

1) 用pw.x进行电子密度的自洽计算;

2) 用ph.x对小的q网格点进行动力学矩阵元的计算; 【得到fildyn='$name.dyn'动力学矩阵文件，在后处理中需要对该文件处理】

以下步骤可以看做是后处理步骤，统一写入phonon_spectra_3.2脚步文件，一次执行

3)用q2r.x计算实空间力常数矩阵；【输入fildyn='$name.dyn'，产生flfrc='$name.fc' (用来设置输出力常数矩阵的文件)；$name.fc被matdyn.x程序调用】

4)用matdyn.x计算出声子色散曲线；【调用flfrc='Sc444.fc'，得到flfrq='$name.freq'(用来输出频率的文件)】

5)用matdyn.x计算声子态密度。【调用flfrc='Sc444.fc', flfrq='Sc444.freq'，得到fldos='$name.dos'：用来存储声子态密度的文件】

以下转载自Valenhou： PWSCF计算晶体的声子色散曲线和态密度

pwscf 是采用线性响应的方法来进行晶格动力学性质的计算。在计算晶体的声子色散和态密度时的步骤：

i)用pw.x进行自洽计算;

ii)用ph.x对小的q点网格进行计算，得到这些q点的动力学矩阵元;

iii)用q2r.x计算出实空间中的力常数矩阵;

iv)用matdyn.x计算声子色散曲线;

v) 用matdyn.x计算声子态密度。

下面以Sc为例子并针对pwscf的最新版本3.2.1来说明（早期版本在计算声子色散曲线较麻烦，因为它不能自动处理q点网格，然后对每个q点一次性计算，而是需要手动产生这些点，一个个计算)。

1) 用pw.x进行电子密度的自洽计算
&control
title='Sc, hexagonal cell'
calculation = 'scf'
restart_mode='from_scratch',
prefix='sc',
pseudo_dir = './',
outdir='./tmp'
tprnfor=.true.
/
&system  
ibrav=4,
celldm(1)=6.05606,
celldm(3)=1.71298,
nat=2,
ntyp=1, nbnd= 30,
ecutwfc=30.0,
occupations ='smearing', degauss =0.01
smearing ='mp'
/
&electrons
diagonalization='cg'
diago_cg_maxiter= 60
mixing_mode = 'plain'
mixing_beta = 0.5
conv_thr = 1.0d-6
/
ATOMIC_SPECIES
Sc 44.955910 Sc.pw91-nsp-van.UPF
ATOMIC_POSITIONS (crystal)
Sc 0.3333333333333286 0.6666666666666714 0.2500000000000000
Sc 0.6666666666666714 0.3333333333333286 0.7500000000000000
K_POINTS (automatic)
8 8 6 0 0 0
注 意Sc是金属，在此例子中，我们选用MP方法来确定电子的占有数(见occupations ='smearing', smearing ='mp')，这里未经测试而选用了展宽系数为0.01 Ry (见degauss=0.01)。在进行声子色散曲线的计算时，不必需对 calculation设置为'phonon'，在新版本中，直接设置为scf。这一步计算产生势以及电荷密度供一下的计算中利用到。【何处用到？】

2)用ph.x对小的q网格点进行动力学矩阵元的计算
phonon for Sc
&inputph
tr2_ph=1.0d-10,
prefix='sc',
fildvscf='scdv',   #File where the the potential variation is written (for later use in electron-phonon calculation).
amass(1)=44.955910,
outdir='./tmp',
fildyn='sc.dyn',
elph=.false.,
trans=.true.,
ldisp=.true.
nq1=4, nq2=4, nq3=2
/
注 意这里trans和ldisp必须设置为.true.。其中trans为.true.表示要计算声子相关的性质，ldisp设置为.true.表示要计算 声子色散曲线。另外 prefix和outdir的设置尽量与上一步自洽计算中的设置一致，以能读入上一步计算得到的数据。另外nq1,nq2和nq3是用来设置q网格点的。 为了得到实空间的力常数矩阵，这里采用的是先计算出q空间中小的q网格点的动力矩阵元，然后采用fft变换得到实空间的力常数矩阵。因此在这一步计算中需 设置小的q网格点的网格大小。【为什么要设置成小的，多小算小？】

3)用q2r.x计算实空间力常数矩阵  【输入name.dyn，产生

name.fc】
&input
zasr='simple', fildyn='sc.dyn', flfrc='Sc444.fc', la2F=.false.
/
在q2r.x的输入文件中需指定
fildyn： 用来设置包含了q网格点的动力学矩阵元的文件，与上一步的 fildyn设置一致;
flfrc：用来设置输出力常数矩阵的文件;
la2F：用来设置是否计算出实空间中电-声耦合系数;针对计算材料的超导性质;
zasr:如何处理‘声学支求和规则“，该规则是用在处理Born有效电荷的，要求Born有效电荷的总和是零。可赋的值有:
no,表示不处理声学支求和问题
simple, 表示通过对力常数矩阵的对角元素进行修正来考虑3支声学横模的求和处理;
crystal,
one-dim,
zero-dim
这里我们设置为simple。

4)用matdyn.x计算出声子色散曲线
&input
asr='simple', amass(1)=44.955910,
flfrc='Sc444.fc', flfrq='Sc444.freq', la2F=.false., dos=.false.
/
131
0.000000 0.000000 0.000000 0.00
0.000000 0.016667 0.000000 0.00
0.000000 0.033333 0.000000 0.00
0.000000 0.050000 0.000000 0.00
0.000000 0.066667 0.000000 0.00
0.000000 0.083333 0.000000 0.00
0.000000 0.100000 0.000000 0.00
0.000000 0.116667 0.000000 0.00
0.000000 0.133333 0.000000 0.00
0.000000 0.150000 0.000000 0.00
0.000000 0.166667 0.000000 0.00
0.000000 0.183333 0.000000 0.00
0.000000 0.200000 0.000000 0.00
0.000000 0.216667 0.000000 0.00
0.000000 0.233333 0.000000 0.00
0.000000 0.250000 0.000000 0.00
0.000000 0.266667 0.000000 0.00
....
这里要输入131个特殊q点的坐标。与计算能带结构时一样，需先选出要计算的高对称q点的走向以及高对称点的坐标，然后产生这些线上的q点的坐标。

计算出来的每一个q点的本征频率可按上一个blog中提到的方法处理一下后画图。

5),用matdyn.x计算声子态密度
&input
asr='simple', amass(1)=44.955910,
flfrc='Sc444.fc', flfrq='Sc444.freq', la2F=.false., dos=.true.
fldos='phonon.dos', nk1=10, nk2=10, nk3=10, ndos=50
/

要计算声子态密度,dos必须设置为.true.另外fldos用来设置输出的态密度值，计算态密度时要用更密的q点网格，这需设置nk1, nk2, nk3。另外还有态密度的能量刻度上的点的数目，由ndos来设置。

手册5.0.2例子1解读

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-PP

cat > 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.UPF

ATOMIC_POSITIONS

Ni 0.00 0.00 0.00

K_POINTS AUTOMATIC

4 4 4 1 1 1

EOF

$ECHO "  running the scf calculation for Ni...\c"

$PW_COMMAND < ni.scf.in > ni.scf.out

check_failure $?

$ECHO " done"

# phonon calculation at X  【这是非自洽计算？】

cat > ni.phX.in << EOF

phonons 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.0

EOF

$ECHO "  running the phonon calculation at X for Ni...\c"

$PH_COMMAND < ni.phX.in > ni.phX.out

check_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   0.0000000), wk =   0.0625000

k(    2) = (   0.5000000 -0.5000000   0.5000000), wk =   0.2500000

k(    3) = (   0.2500000  0.2500000   0.2500000), wk =   0.2500000

k(    4) = (   0.7500000 -0.2500000   0.7500000), wk =   0.7500000

k(    5) = (   0.0000000 -1.0000000   0.0000000), wk =   0.1875000

k(    6) = (   0.0000000  0.0000000   0.5000000), wk =   0.1250000

k(    7) = (   0.5000000 -0.5000000   1.0000000), wk =   0.2500000

k(    8) = (   0.0000000 -1.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文件中，按格式修改相关参数，执行程序，得到….  运行read_phonon.x，得到可用于绘图的.dat文件【注意采用该程序提取的频率的单位为cm-1，不是THz】

1. pw.x 和ph.x脚本

#!/bin/bash

####################################################################

./home/xiaoqiugood/pwscf/pw/environment/environment_1

exportPARA_PREFIX='mpirun  -np 12 '

# or   export PARA_PREFIX='mpirun' ,export  PARA_POSTFIX= -np 3

exportTMP_DIR=/home/xiaoqiugood/pwscf/pw/tmp/tmp3

rm -r $TMP_DIR/*.save

rm $TMP_DIR/*

export name='a225'

for a in 1000

do

cat >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=2,

celldm(1)=7.836769829,

nat= 4, ntyp= 2

ecutwfc =60,

nbnd= 30,

occupations ='smearing',

smearing ='mp',

degauss =0.01

/

&electrons

mixing_beta = 0.7

conv_thr = 1.0d-8

/

ATOMIC_SPECIES

H     1.008     H.pbe-van_ak.UPF

Sc    44.9559   Sc.pbe-nsp-van.UPF

ATOMIC_POSITIONS(crystal)  #必须用pwscf优化后的原子位置，不能直接用VASP优化后的原子位置

H  0.750000000000.75000000000  0.75000000000

H  0.250000000000.25000000000  0.25000000000

H 0.500000000000.50000000000  0.50000000000

Sc 0.00000000000 0.00000000000 0.00000000000

K_POINTS {automatic}

24 24 24  0 0 0

EOF

PARAPREFIX

PW_ROOT/pw.x  <name.scf.in_a >name.scf.out_

a

done

# occupations用来设置确定电子占有数的方法，赋值为'smearing'表示采用smearing的方法来确定电子的占有数，随后须设置smearing和degauss关键词。 smearing用来指明确定电子占有数的一种具体的smearing方法，赋值为'gaussian'表示采用Gaussian函数来确定电子占有数；赋值 'mp'，即表示采用 Methfessel-Paxton first-orderspreading (see PRB 40, 3616 (1989))方法或函数来确定电子占有数。

degauss用来确定smearing方法中有关函数的展宽参数，赋值为0.05表示设置上面函数中的展宽参数为0.05。degauss 为0，相当于fixed-occupation的计算，即每个态的电子占有数是固定的。当体系为半导体或绝缘体时，可以设置degauss=0，其他情况下不能。

# Variable: nbnd;Type:  INTEGER  Default:        for an insulator, nbnd = number of valence bands (nbnd=nelec/2,see below for nelec); for a metal, 20% more (minimum 4 more)  Description: number of electronic states(bands) to be calculated. Note that in spin-polarized calculations the numberof k-point, not the number of bands per k-point, is doubled

cat>$name.ph.in << EOF

phonon of$name  at Gamma

&inputph

tr2_ph=1.0d-12,

prefix='$name',

amass(1)=1.008

amass(2)=44.955910,

trans=.true.,

ldisp=.true.,

outdir='$TMP_DIR',

fildyn='$name.dyn',

nq1=4,nq2=4,nq3=4 ##该参数如何设置？   不只是gamma点，那是哪些点，高对称性点包含了吗？

# 本部分重点设置参数：tr2_ph（频率收敛截断）, fildyn动力学矩阵文件）, amass（原子质量）。如果是绝缘体还要考虑加入玻恩有效电荷的计算参数：epsil=true。

# tr2_ph . Threshold forselfconsistency.Default: 1e-10.可以调高些，计算量增加的不多，但是能够提高声子的计算精确度。

#注意这里trans和ldisp必须设置为.true.。其中trans为.true.表示要计算声子相关的性质，ldisp设置为.true.表示要计算声子色散曲线。

手册解释：Description: if .true. the phonons are computed. if trans .and. epsil effective charges are calculated.  Default: .true.

手册解释ldisp：Default:  .false. Description:   If .TRUE. the runcalculates phonons for a grid of  q-points specified by nq1, nq2, nq3  - for direct calculation of the entire phonon dispersion.  The pw.x data file should not be producedusing "calculation='phonon'" in this case.

# fildyn='$name.dyn'：这个是输出的动力学矩阵文件，在后处理中需要对该文件处理，大家可以看看里面是什么内容。手册介绍：File where the dynamical matrix is written. Default:        'matdyn'.

# 另外 prefix和outdir的设置尽量与上一步自洽计算中的设置一致，以能读入上一步计算得到的数据【哪一个参数是读入上一步数据】。

#nq1,nq2和nq3是用来设置q网格点的。为了得到实空间的力常数矩阵，这里采用的是先计算出q空间中小的q网格点的动力矩阵元，然后采用fft变换得到实空间的力常数矩阵。因此在这一步计算中需设置小的q网格点的网格大小。【q点选择依据，怎么不见设置高对称性q点，高对称性q点与k点是一致的吗？<20】手册解释：Variables:  nq1, nq2, nq3. Default: 0.  Description: Parameters of the Monkhorst-Pack grid (no offset) used when ldisp=.true.Same meaning as for nk1, nk2, nk3 inthe input of pw.x.  【用ldisp=.true.自动产生q点?】

# 注意Variables: xq(1)  xq(2)  xq(3)设置在inputph namelist之外；而iq1, iq2, iq3 或nq1、nq2、nq3均设置在namelist之内。These (iq1, iq2, iq3)go together with nq1, nq2, nq3 and allow to choose just one point out of theMonkhorst-Pack grid with ldisp=.true. Note the the actual point chosen is something like  (iq1-1)/nq1, (iq2-1)/nq2, (iq3-1)/nq3 (so, check the outputfor what you get). Also make sure that PW left *.wfc)。

#手册解释Variables:xq(1)  xq(2)  xq(3). Decription:    The phonon wavevector; must be equal to the one usedin the non-selfconsistent calculation (not read if ldisp is true).

# 关于高对称性点，与计算能带类似，要自己先选定一些高对称点，并产生这些高对称点之间其他点。

/

EOF

PARAPREFIX

PW_ROOT/ph.x < name.ph.in>

name.ph.out

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.0000000E+00           1

-0.9833333      0.5000000      1.6666669E-02           2

-0.9666667      0.5000000      3.3333331E-02           3

.......(略)

-0.7666667      0.7333333      0.0000000E+00         149

-0.7583333      0.7416667      0.0000000E+00         150

-0.7500000      0.7500000      0.0000000E+00         151

EOF

$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：用来存储声子态密度的文件;

# deltaE表示energy grid step (eV)

例子文件：

# job script for quantum espresso

# the name of queue, should be cscore or csnode

APP_NAME="cscore"

# number of processes

NP=16

# fix command style

RUN="RAW"

# current directory

CURDIR=$PWD

# build .nodelist

rm -rf $CURDIR/.nodelist >& /dev/null

for i in `echo $LSB_HOSTS`

do

echo i>>

CURDIR/.nodelist

done

# directory for quantum espresso

EXEDIR=$HOME/bin/espresso502/bin

# directory for mpi (mvapich)

MPIDIR=/home/compiler/mpi/mvapich/1.0/icc.ifort-11.1/bin

# start calculation (scf)

#MPIDIR/mpirunnp

NP -machinefile CURDIR/.nodelist

EXEDIR/pw.x -in scf.in > scf.out

# start calculation (nscf)

#MPIDIR/mpirunnp

NP -machinefile CURDIR/.nodelist

EXEDIR/pw.x -in nscf.in > nscf.out

# start calculation (phonon)

#MPIDIR/mpirunnp

NP -machinefile CURDIR/.nodelist

EXEDIR/pw.x -in ph.in > ph.out

#PARA_PREFIX='MPIDIR/mpirunnp

NP -machinefile $CURDIR/.nodelist'

#export PARA_PREFIX

####################################################################

PW_ROOT=/home/users/spclibtt/bin/espresso502/bin

PSEUDO_DIR=/home/users/spclibtt/home/yexq/shanghai-opt/pw-work/pseudo

TMP_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 3

export PATH=/home/compiler/mpi/mvapich/1.0/icc.ifort-11.1/bin:$PATH

export PARA_PREFIX='mpirun -np 16'

export PW_ROOT PSEUDO_DIR TMP_DIR

export name='sch4-i4mm'

rm -r $TMP_DIR/*.save

rm $TMP_DIR/*

###########self-consistent calculation ######################

for a in 2500

do

cat > 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_SPECIES

H     1.008      H.pbe-van_ak.UPF

Sc    44.9559    Sc.pbe-nsp-van.UPF

CELL_PARAMETERS {bohr}

4.745568648     -4.27115E-08    0

-4.27115E-08    4.745568539     -9.49146E-09

4.74573E-09     -1.89829E-08    8.842464528

ATOMIC_POSITIONS {crystal}

H        0.499999965   0.500000040   0.370850743

H        0.499999950   0.500000054   0.629149579

H       -0.000000022   0.000000030   0.870851011

H       -0.000000027   0.000000036   0.129148779

H       -0.000000013   0.500000019   0.249999381

H       -0.000000012   0.500000019   0.750000700

H        0.499999968   0.000000027   0.249999365

H        0.499999976   0.000000024   0.750000717

Sc      -0.000000006  -0.000000010   0.499999832

Sc       0.500000021   0.500000003   0.000000162

K_POINTS {automatic}

16 16 9 0 0 0

EOF

MPIDIR/mpirunnp

NP -machinefile CURDIR/.nodelistEXEDIR/pw.x  < name.scf.in_a >name.scf.out_

a

done

cat >$name.ph.in << EOF

phonon of $name  at  X

&inputph

tr2_ph=1.0d-12,

prefix='$name',

fildvscf='$name.dv',

amass(1)=1.008

amass(2)=44.955910,

trans=.true.,

ldisp=.true.,

outdir='$TMP_DIR',

fildyn='$name.dyn',

nq1=4,nq2=4,nq3=2

/

EOF

MPIDIR/mpirunnp

NP -machinefile CURDIR/.nodelistEXEDIR/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 3

export  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.out

cat > 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.0000000E+00           1

-0.9833333      0.5000000      1.6666669E-02           2

-0.9666667      0.5000000      3.3333331E-02           3

.......

-0.7500000      0.7500000      0.0000000E+00         151

EOF

$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

原作者：叶小球

http://blog.sciencenet.cn/home.php?mod=space&uid=567091