5.0.2版计算电声耦合常数分为以下几个步骤：

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 el-phon coefficients (inparticular the double-delta  integral at Ef is very critical).

问题1：第一步就涉及q点和k点设置问题，例子中是设置成4 ：2 ：1的关系，即：16x16x16致密k点自洽计算, 8x8x8, 4x4x4。请问q点的设置依据是什么？一般应设置成多少合适？
       什么是 thedouble-delta  integral at Ef ？ 

   If nq1,nq2,nq3 are not divisors of that numbers used to generate dense k-points set, 则会报错。

查看elph.f90源程序：

call cryst_to_cart (1, xq, at, -1)

 qx = nint(nk1fit*xq(1))   #NINT是四舍五入取整

 if (abs(qx-nk1fit*xq(1)) > eps) &

      call errore('elphsum','q is not a vector in the dense grid',1)

2. make aself-consistent calculation for Alusing a grid of k-points   that is suitable for good self-consistencyand phonon calculation.This example uses a (888) MP grid.


3.  make the phonon andelectron-phonon calculation for the grid of q-points.   Specifyelph=.true., and the name of a file where the derivative of the  potential is stored "fildvscf".This example uses a (444) MPgrid of   q-points (nq1=4, nq2=4, nq3=4).
  The output contains the resultsfor the el-ph coefficient at each q-point  lambda(q), gamma(q), andthe double-delta integral at several values of   the gaussianbroadening (set in filePH/elphon.f90). These are useful   forconvergence testing.
 (input=al.elph.in,output=al.elph.out)
  The results are also written intooutputfiles "a2Fq2r.*", one per value   of the gaussianbroadening, for further use.
 
问题2：第三步中（1）5.0.2新版在ph.x步骤设置 electron_phonon='interpolated',  以前的老版本是设成：elph=.true.,请问两个参数的区别在哪里？

    （2）convergence testing是算完之后测试？  gaussian broadening 对计算结果有什么影响？

4. Bring to r-space both force constants andel-phon coefficients using "q2r.x"   Output in files"a2Fmatdyn.*", one per value of the gaussian broadening
问题3：第四步q2r.x，为什么要转换到实空间？与gaussian broadening 有什么关系？这一步是必须的吗？ 

5.  Calculate gamma  onselected lines using "matdyn.x" (dos=.false.)
问题4：这一步是在计算声子色散曲线，它是计算电声耦合常数必须的步骤吗？

6.  Calculate  lambdacoefficient (in file"lambda" and the a^2F(omega)
  function using"matdyn.x" (dos=.true.)

问题5：这一步感觉是在计算声子态密度，请问该步能同时给出 lambda coefficient和a^2F(omega)
  function 吗？
7.  Calculate lambda coefficient (in"lambda.out"and Tc using"lambda.x"

问题6：计算仍在进行中，请问有过计算经历的同仁，新版通过这一步可以直接算出Tc吗？

       综上，感觉第5步，第6步似乎不是计算电声耦合常数必须的步骤，不知理解是否正确。完整的例子文件略。

电声问答：

 Q::

是的，是因为这个问题，后来楼上的师兄指导把elph文件里面的负值改成0就可以产生了Tc了，不会出现NaN,不过有意思的是，之所有有两个elph里面是负值居然是因为计算的ph.in里面的精度不够，之前的精度是10-12，现在是-14这下就没有负值了， 不过这计算的时间也是相当伤不起，到今天也没有算完，中途还断了一次电，唉。
对了，请问，您最后看lambda.out结果时是看差不多收敛就可以了，还是要测试degauss，因为最后结果会从degauss=0.005-0.05，暑假在北京学习时一位师兄说算之前选一个q点计算进行测试，然后直接选测试结果的degauss对应的lambda和Tc，谢谢！

A:::

the smallest（broadening） that yields converged values for the number of k-points you can afford.

其实保罗这句话挺含糊的，我觉得是取文件中这些值里面让lambda收敛的最小的展宽，并且你这个展宽在你当前k点所允许的范围内。如果展宽太小，k点不密的话，在整个能量空间中会出现很多尖峰，导致计算错误。所以我猜他是想说这个。

这个文献里面的值也是这个值，结论应该和您的结论差不多，文献是：PRL     101, 026403 (2008) 师兄的测试的方法我还在算，还没有得到结果，师兄的方法是这样的：我算一个密集K点自洽和一般K点的自洽，然后就单独算一个q点的电声，当然这个q点不能是gamma点，在计算后会得到elph.0.0000 0.50000 0.0000类似的这样计算的q点的文件，里面会对应每个broadening（0.005-0.05）下面都有好多lambda值，把每个broadening下面的lambd值进行相加，就可以得到一个横坐标是broadening，纵坐标是lambda的一条关系曲线，然后再通过改变自洽里面的degauss，就可以得到每一个degauss对应一条这样的关系曲线，测试应该会有趋于收敛的，到时候就可以取那个degauss进行整体的计算，然后最后取那个broadening对应的lambd和Tc。

我讲计算结果处理了下，其中442一共12个q点，然后选取第11个q点作为测试q点，因为它计算最快，然后给出测试结果，如图。如果按照图中的话，那么在自洽的时候可以degauss可以为0.01或者0.02，在lambda.out中，结果应该选择0.03所对应的结果，也就是lambda=0.16。对了，上面的回复中，您提到，lambda过小，电子和声子不太喜欢耦合在一起，我看了下example里面的Al的lambda也差不多0.42左右，这个lambda值的大小有个什么范围么，就是多少范围属于大的，可以耦合，那些范围是太少，不太喜欢耦合？谢谢！

test.jpg

A::::

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

Q:

Gaussian Broadening:   0.050 Ry, ngauss=   0
    DOS = 44.226975 states/spin/Ry/Unit Cell at Ef= 14.581045 eV
    lambda( 1)=  0.0057   gamma=    0.37 GHz
    lambda( 2)=  0.0057   gamma=    0.37 GHz
    lambda( 3)=  0.0057   gamma=    0.37 GHz

....................

A:::

gamma的频率是电声耦合导致的phonon linewidths，貌似你自己已经回答了

A:::

1、在lambda.x的input文件的第一行，需要输入3个参量，分别是emax，degauss，以及smearing method。请问这三个参数有什么含义，取法有什么讲究，对于计算的结果有什么影响？

2、在lambda.x的output文件里面，分别列出了degauss从0.005到0.05的lambda及Tc的计算结果，不同的degauss对应的结果不同，那么，到底哪一个结果才是我们想要的最终结果呢？

Q:::

emax是比最高声子能量还高一点的一个值，取多少自己定。声子能量，直接看最高声子频率就是，要取以THz为单位的。
degauss和smearing method是在布里渊区内取点的时候用的展宽方式和大小，但具体该怎么取我还不懂
因为你布里渊区取点，k点数目是一定的，那么degauss越大越容易收敛，所以你就看degauss比较大的几个lambda和Tc，是否随着degauss的增加有收敛的趋势。如果有这种趋势，那么恭喜你。

，最后一行取0.1到0.15之间，是代表电子之间的库仑排斥作用

A::::

Dear Elie,>from elphsum : error #         2>     q is not a vector in the dense gridThat means k-points set used in the first scf step and what you specified in elph.in file are incompatible, i.e. nq1,nq2,nq3 are not divisors of that numbers used to generate dense k-points set.Regarding lambda.x, I guess, there was an example file (an input file) where an explanation was done.  If you can not find it (please search before) then contact me.Besides, all info calculated via lambda.x  can now be obtained using matdyn.xBests,Eyvaz.

A:::: 16:48:44

      PWSCF 5.0.2版计算电声耦合常数时在ph.x步骤设置 electron_phonon='interpolated',   以前的老版本是设成：elph=.true.,请问两个参数的区别在哪里？

Q::::::::::::48776252) 16:49:39

如果你设置electron_phonon＝.true. 两者是一样的。

QQQ::::

    新版本的意思， 是用稀疏的k点算phonon

    dense的k点算elph coupling

A:::::::

    是用稀疏的k点算phonon, dense的k点算elph coupling，与下面那个步骤对应：

1）选取很密的k点网格，用pw.x进行自洽计算；

2）、选取略疏的k点网格，用pw.x进行自洽计算

3)、选取较疏的q点网格，用ph.x进行声子的计算，得到每一个q点的动力学矩阵；

如：Electron-phonon coefficients for Al

&inputph

tr2_ph=1.0d-10,

prefix='al',

fildvscf='aldv',

amass(1)=26.98,

outdir='./tmp' ,

fildyn='al.dyn',

elph=.true.,   #5.0.2改为electron_phonon='interpolated',

trans=.true.,

ldisp=.true.

nq1=4, nq2=4, nq3=4

/

A:::::     第三步中q点的数量一般不少于50个？

Q:::::

其实电声还是只适合高对称性结构， 低对称性的很费劲。

tr2_ph=1.0d-12  

一般是k点的一半就可以。

Q:::::

  就是，感觉电子和声子要发生相互作用的话，声子得有很多才行，即q点应该取很多

Q::::: 恩，当然。

A:::::

     谢谢。另外，感觉后面的步骤都应该属于后处理的步骤，而且第5、6步计算声子色散曲线和态密度，对于计算电声耦合常数似乎是不必要的步骤：

4)、 用q2r.x由每一个q点的动力学矩阵得到实空间中力常数矩阵Bring tor-space both force constants and el-phon coefficients using "q2r.x"

5)、声子色散曲线计算Calculate gamma  on selected lines using "matdyn.x" (dos=.false.)

6)、声子态密度计算Calculate  lambda coefficient (in file"lambda") and the a^2F(omega)  function using "matdyn.x" (dos=.true.)  

7)   Calculate lambda coefficient (in "lambda.out")and Tc using "lambda.x"

Q:::::

yes

A:::::

        Thanks. 自洽计算步骤中degauss =0.2参数的设置对电声耦合计算结果有什么影响，需要测试吗？

Q:::

是scf吗？ 0.2太大了。  一般0.05－0.03就可以

附--大师10年前的提问：

[Pw_forum] Question for electron phonon coulpingStefano de Gironcoli 
Tue Sep 23 11:41:14 CEST 2003

• Previous message: [Pw_forum] Question for electron phonon coulping

• Next message: [Pw_forum] compilation error

• Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

Yanming mawrote:>>Do a calculation with a given set of k points and calculateel-ph as a>>function of sigma>>>In practise,>I choose a 8x8x8 k mesh and a q point (say 0.5 0.5 0.0). I set nsig >=4, 8, 12, 16, 20, 24...respectively. I calculate the lambdas at this >q point with this k mesh and with different nsigs. Concerning the >lambda value, For example, for nsig=4 and 8, the final lambda values >are choosen as the final step gaussian broadening (step 4 for nsig=4, >step 8 for nsig=8).>Then I can get a curve for nsig.vs.lambda at this q point.You do not need to do a separate calculation for each value of nsig. with the original value of nsig=10 and the  step size of 0.01 (ryd) in a sigle run you get results for degauss1=0.01, 0.02, 0.03,,..., 0.10 (ryd).The result of step 4 for nsig=4 is identical to step 4 for nsig=8,12,16,20,24 or any N>4 !If you increase nsig to 24 (leaving the step size equal to 0.01) you can reach degauss1=0.24 if you like but, remember, one is interested in the limit of degauss1 going to zero.NOTICE that I call degauss1 the value of the smearing used in the calculation of el-ph interaction which  will differ in general from degauss used in the scf/phonon calculation.>>Repeat the calculation with a larger k mesh, for example, 12x12x12.>>I compare the two curves and get a reliable nsig value, which the two >curves start to agree.>>if I repeat the calculation with more denser k mesh (say 16x16x16), I >will get a smaller reliable nsig by comparing the two more denser k mesh.>>>Trust only the curve in the region in which the calculation is converged>>with respect to k-ponts summation.>>In practise, how can I decide the calculation is converged with >respect to k points summation?>At a q point, with the same nsig, The calculated lambdas are same >within some errors for two K points mesh. Does this mean convergence? When for a given value of degauss1 (a given step in the nsig loop) the calculated lambda does not change anymore with increasing density of the  K-points mesh then you have reached convergence for that value of degauss1.>>>For sufficiently large set of k-points the SAFE region will approach>>enough the limit of sigma going to zero>>that you can stop the procedure.>>>>I can not understand this paragraph. Sigma going to zero...?>SAFE meant converged with respect to k-point mesh.The idea is that for sufficiently large set of K-points even the smallest values of degauss1are converged and one can use these values as approximations of the degauss1 going to zero limit.How close one needs to go to the limit depends on the shape of the lambda vs degauss1curve.Let's say 0.01-0.03 ryd or larger if the curve is sufficiently flat (horizontal).As an example of this idea in a different context I include a figure where convergence with respect to k-points and degauss is computed for the total energy of lead.In this figure one sees that for ngauss=1 the "converged" courve is sufficiently flat that you can usea value of degauss as large as 0.05 ryd being within a franction of mryd from the exact limit (degauss=0 and infinite number of k-points). ngauss=0 converge to the same limit but with a stronger dependence on degauss.This is the reason why in the calculation of lambda ngauss1 is always chosen equal to 1, although for el-ph interaction the dependence on degauss1 might be less nice.So, do the calculation for several  set of k points and a given value of nsig (24 if you like). This will give you an entire curve lambda(degauss1) for each set of k-points.Increse the k-point grid until the region of converged degauss1 extent sufficiently close to zero that you feel safe in extrapolating  that value to zero.How much you need to go will depend on the shape of the function: whether it is very flat as the ngauss=1 case in the figure or with a large and changing slope as the negauss=0 case.Once you are satisfied with the convergence of lambda at a given value of phonon wavevector q, then you still need to check convergence with respect to number of q-points (phonon wavevectors) you include in the calculation.So you should compute el-ph interaction for a larger set of q-points and REPEAT the above procedure because theconvergence properties will in general depend on the q-points....... and so on an so forth ...Stefano de Gironcoli-------------- next part --------------A non-text attachment was scrubbed...Name: EN_VS_SIGMA.psType: application/postscriptSize: 7073 bytesDesc: not availableUrl : /pipermail/attachments/20030923/68af8f37/attachment.ps

//////////////////////////////前一篇提问////////////////////////////////////

[Pw_forum] Question for electron phonon coulpingYanming ma 
Mon Sep 22 20:18:21 CEST 2003

• Previous message: [Pw_forum] Re: Ga 3d core electrons as valence electrons (fhi2upf)

• Next message: [Pw_forum] Question for electron phonon coulping

• Messages sorted by: [ date ] [ thread ] [ subject ] [ author ]

Dear Stefano de Gironcoli,Thanks very much for your excellent comments.I am not sure if I can understand 100 percent of your comments .I try to follow your comments in order to make sure I can do the calculation in a right way.>Do a calculation with a given set of k points and calculate el-ph as a>function of sigmaIn practise,I choose a 8x8x8 k mesh and a q point (say 0.5 0.5 0.0). I set nsig =4, 8, 12, 16, 20, 24...respectively. I calculate the lambdas at this q point with this k mesh and with different nsigs. Concerning the lambda value, For example, for nsig=4 and 8, the final lambda values are choosen as the final step gaussian broadening (step 4 for nsig=4, step 8 for nsig=8).Then I can get a curve for nsig.vs.lambda at this q point.>Do again the same calculation increasing the number of k-points and the>same set of values for sigma.Repeat the calculation with a larger k mesh, for example, 12x12x12.>You will notice that for sufficiently large values of sigma the two>curves will agrees. This means that>the k-points mesh is sufficiently accurate for sigmas larger than that>but may not be good for smaller values.I compare the two curves and get a reliable nsig value, which the two curves start to agree.>Repeating the calculation for even larger set of k-points (same set of>sigma) you will see that from a value>of sigma that is smaller than the previous one the most accurate>calculation (the one using more k-points)>will agree with the second-best (and with the first, less accurate,>calculation when these last two agree)if I repeat the calculation with more denser k mesh (say 16x16x16), I will get a smaller reliable nsig by comparing the two more denser k mesh.>Trust only the curve in the region in which the calculation is converged>with respect to k-ponts summation.In practise, how can I decide the calculation is converged with respect to k points summation?At a q point, with the same nsig, The calculated lambdas are same within some errors for two K points mesh. Does this mean convergence?>For sufficiently large set of k-points the SAFE region will approach>enough the limit of sigma going to zero>that you can stop the procedure.>I can not understand this paragraph. Sigma going to zero...?Any comments?Thanks in advance.Yanming Ma

事例：

####################1st a scf calculation for the dense k-point  grid,specify   option la2f=.true #############################

cat > $name.scf.fit.in << 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= 70, ecutrho = 600.0

 occupations ='smearing', smearing ='mp',degauss =0.02

 la2F = .true.

/

&electrons

   mixing_beta = 0.7

   conv_thr =  1.0d-7

/

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}

20 20 11 0 0 0

EOF

$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/pw.x  < $name.scf.fit.in >$name.scf.fit.out

############################2nd: a normal scf + phonon dispersion calculation on the coarse k-point   grid, specifying option electron_paononh='interpolated'. and the file name : fildvscf###################

cat > $name.scf.in << 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= 70, ecutrho = 600.0

 occupations ='smearing', smearing ='mp',degauss =0.02

/

&electrons

   mixing_beta = 0.7

   conv_thr =  1.0d-7

/

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/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/pw.x  < $name.scf.in >$name.scf.out

######################################phonon calculation, specifying option elph=.true.#########################

cat > $name.elph.in << EOF

Electron-phonon coefficients for ScH4

&inputph

tr2_ph=1.0d-12,

prefix='$name',

amass(1)=1.008

amass(2)=44.955910,

fildyn='$name.dyn',

fildvscf='$name.dv',

electron_phonon='interpolated',

trans=.true.,

ldisp=.true.,

nq1=4,nq2=4,nq3=2

/

EOF

$MPIDIR/mpirun -np $NP -machinefile $CURDIR/.nodelist $EXEDIR/ph.x  < $name.elph.in > $name.elph.out